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ABSTRACT 


Energy release rate is usually suggested as a quantifier for assessing structural 
damage tolerance. Computational prediction of energy release rate is based on composite 
mechanics with micro-stress level damage assessment, finite element structural analysis 
and damage progression tracking modules. This report examines several issues associated 
with energy release rates in composite structures as follows: Chapter I demonstrates 
computational simulation of an adhesively bonded composite joint and validates the 
computed energy release rates by comparison with acoustic emission signals in the 
overall sense. Chapter II investigates the effect of crack plane orientation with respect to 
fiber direction on the energy release rates. Chapter III quantifies the effects of 
contiguous constraint plies on the residual stiffness of a 90° ply subjected to transverse 
tensile fractures. Chapter IV compares ICAN and ICAN/JAVA solutions of composites. 
Chapter V examines the effects of composite structural geometry and boundary 
conditions on damage progression characteristics. 
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Chapter I. Energy Release of Adhesively Bonded Joints 


In this chapter, the adhesive bond strength of lap-jointed graphite/aluminum composites is 
examined by computational simulation. Computed micro-stress level energy release rates are used to 
identify the damage mechanisms associated with the corresponding acoustic emission (AE) signals. 
Computed damage regions are similarly correlated with ultrasonically scanned damage regions. 
Results show that computational simulation can be used with suitable NDE methods for credible in- 
service monitoring of composites. 

KEY WORDS: Computational simulation, Energy release rate. Nondestructive Evaluation (NDE) 

1.1 INTRODUCTION 

For effective structural health monitoring, it is important to quantify damage tolerance of a 
candidate structure. Since continuous fiber composites are able to arrest cracks and prevent self-similar 
crack propagation, composite structures have received a great deal of consideration for design with 
emphasis on damage tolerance. However, a number of design parameters such as fiber orientation 
patterns, choices of constituent material combinations, ply drops and hybridization, result in complex 
design options for composite structures. Thus, it is necessary to evaluate damage initiation in a 
composite structure and its fracture propagation characteristics for achieving a rational damage tolerant 
design. 


Compared with homogeneous materials, damage initiation and progression characteristics of 
fiber composites are much more complicated. Composite structures often contain some pre-existing or 
induced flaws in matrix and fibers after fabrication of composites. At lower stresses, matrix is likely to 
be cracked because of flaw-induced stress concentrations and cause the matrix flaws to propagate 
across the composite. With the use of established material modeling and finite element models, and 
considering the influence of local defects, through-the-thickness cracks and residual stresses, 
computational simulations have made it possible to evaluate the details of progressive damage and 
fracture in composite structures. In a computational simulation, damage evolution quantifier such as 
the damage volume, exhausted damage energy, and the damage energy release rate (DERR) are used 
to quantify the structural damage tolerance at different stages of degradation. Low DERR levels 
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usually indicate that degradation takes place with minor resistance by the structure. Structural 
resistance to damage propagation is often dependent on structural geometry and boundary conditions 
as well as the applied loading and the state of stress. 

In certain cases such as the room temperature behavior of composites designed for high- 
temperature applications, internal damage initiated as microcracks in the matrix become enlarged to be 
externally visible. Thus, matrix cracking and its effect on damage propagation/damage tolerance need 
be evaluated. Some simulations [1,2] have been successful in predicting damage tolerance and failure 
load of composite structure by considering ply stresses and the corresponding stress limits for matrix 
crack growth. In this report, lap-jointed composite specimens subjected to uniaxial tension are 
investigated. Damage initiation, growth, accumulation, and propagation to fracture are studied. Since 
the complete evaluation of ply and subply level damage/ fracture processes is the fundamental premise 
of computational simulation, a microstress level damage index is added for the identification and 
tracking of subply level damage processes. Computed damage regions are similarly correlated with 
ultrasonically scanned damage regions. Simulation is validated by comparison with test data from 
acoustic ultrasonic (AU) testing [7], Results show that computational simulation can be used with 
suitable NDE methods for credible in-service monitoring of composites. 

1.2 METHODOLOGY 

Computational simulation is implemented via integrating three modules: (1) composite 
mechanics, (2) finite element analysis, and (3) damage progression tracking. The composite mechanics 
module (Murthy and Chamis 1986) is designed to analyze fiber composite structures with an updated 
composite mechanics theory. Its main function is to calculate ply and composite properties of laminates 
from the fiber and matrix constituent characteristics and the composite lay-up. Additionally, it 
determines the composite structural response and ply stresses from the FEM analysis results. In 
simulation, the composite mechanics module is called before and after each finite element analysis. 

The finite element analysis module is able to process linear and nonlinear static and dynamic 
analysis. Four— node anisotropic thick shell elements are usually used to model laminated composites 
(Nakazawa et al 1987). The finite element analysis module accepts laminate properties from the 
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composite mechanics module and performs the structural analysis at each load increment. After 
structural analysis, the computed generalized node stress resultants and deformations are provided to 
the composite mechanics module. The composite mechanics module computes the developed ply 
stresses for each ply and checks for ply failure modes at each node. Failure criteria applied to detect 
ply failures are based on the maximum stress and modified distortion energy (MDE) criteria for 
combined stress effects (Murthy and Chamis 1986). 

The overall evaluation of composite structural durability is carried out in the damage 
progression module (Minnetyan et al 1990) that keeps track of composite degradation for the entire 
structure. The damage progression module relies on the composite mechanics module for composite 
micromechanics, macromechanics and laminate analysis, and calls the finite element analysis module 
for global structural analysis. If excessive damage is detected, the incremental loads are reduced and 
the analysis is restarted from the previous equilibrium stage. Otherwise, if the increment of loads is 
acceptable, another finite element analysis is performed but the constitutive properties and the finite 
element mesh are updated to account for the damage and deformations from the last simulation. 
Simulation is stopped when global structural fracture is predicted. 

1.3 METHOD OF SIMULATION 

The matrix in orthotropic composite plies is divided into two parts: regions A and B. Region A 
represents the area in which stress concentrations induced by the interaction of matrix and fiber do not 
create any effect in matrix. Region B represents the interaction zone between fiber and matrix. Figure 
1.1 shows the details of regional subdivision in transverse and normal directions of a composite ply 
with square packing. 

Considering the behavior of longitudinal stress cr ni , transverse stress cr l22 , in-plane shear stress 
cr ;i2 , out-of plane shear stress cr /23 , temperature gradient AT,, and moisture M ; , Murthy and Chamis 
(1986) present the complete set of equations for evaluating ply microstresses in regions A and B. For 
example, ply microstresses due to <j iu are given by: 
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transverse stress in region A, and is the matrix transverse stress in region B. If the ply is subject 

to combined stresses, its microstresses are obtained by simply superimposing results of all 
corresponding stress components. Ply transverse fractures usually begin in region B due to the elevated 
stress levels from stress concentration. Microstress level damage tracking is able to quantify the type of 
damage in the matrix by comparison of microstresses with constituent stress limits. A microstress 
damage index is defined as a binary number with 14 bits in the damage progression module. The 14-bit 
number corresponds to the following components: (1) SM1A(+) longitudinal stress in region A, 
positive; (2) SMIA(-) longitudinal stress in region A, negative; (3) SM2A(+) transverse stress in region 
A, positive; (4) SM2A(-) transverse stress in region A, negative; (5) SM2B(+) transverse stress in 
region B, positive; (6) SM2B(-) transverse stress in region B, negative; (7) SM3A(+) normal stress in 
region A, positive; (8) SM3A(-) normal stress in region A, negative; (9) SM3B(+) normal stress in 
region B, positive; (10) SM3B(-) normal stress in region B, negative; (11) SM12A in-plane shear stress 
in region A; (12) SM12B in-plane shear stress in region B; (13) SMI 3 A out-of-plane shear stress in 
region A; (14) SM23A out-of-plane shear stress in region A. When the binary bit corresponding to 
SM2B(+) is set equal to 1, it indicates that region B is fractured with the transverse failure mode. In a 
subsequent stage transverse cracking will spread to region A. 


Mital (1993) presents a more refined slice-by-slice substructuring model for assessing ply 
failure modes at the microstress level. In his model, fiber is substructured into several slices and the 
micromechanics equations are applied at the slice level. Once the equivalent slice properties are 
defined, the ply properties are obtained using classical lamination theory. By using maximum strength 
and combined stress criteria on the fiber, matrix, and interface microstresses, failure criteria are applied 
directly at the slice microstress level. With this model, the variation of interfacial bonding around the 
fiber circumference can be depicted more clearly. As a result, the local matrix cracking and fiber 
breaks can be monitored more closely in computational simulation. For ceramic matrix composites 
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structural life prediction is quantifiable more accurately with the slice-by-slice ply substructuring 
model. The fiber substructuring and slice geometry are shown in Figure 1 .2. 

1.4 SIMULATION RESULTS AND DISCUSSION 

In this report, two unidirectional graphite/aluminum composite plates with single lap joint 
under uniaxial tension are used to demonstrate the use of microstress damage index in computational 
simulation. The specimen has a length of L = 76.2mm (3 in.), a width of W = 19.05mm (0.75 in.) and 
a thickness of H = 1.27mm (0.05 in.). The fiber volume ratio is V f = 0.60 and the void volume ratio is 

V v = 0.01 . The specimens are bonded with adhesive to produce a bond area of 0.75"x0.75". The first 

adhesive type is an epoxy resin and the other is a graphite/epoxy prepreg tape. The specimen is 
assumed to be dry with zero moisture content. Additionally, two metal plates with dimensions of 
2"x0.75" are also bonded to the specimens with epoxy resin. The fiber and matrix properties used for 
computational simulation are given in Table 1.1 and Table 1.2, and the configuration for adhesive 
bonded specimen is shown in Figure 1.3. 

The finite element model as shown in Figure 1.4 has 451 nodes and 360 Mindlin type 
rectangular elements. To represent the test setup described by Quattlebaum (1997), nodes along the 
end of the specimen are restrained against translation to model the fixed boundary and nodes at the end 
of another specimen are constrained only moveable in the longitudinal direction. Axial tension load is 
applied at the moveable end. Moreover, nodes along the moveable end are tied by duplicate node 
constraints to enforce the uniform displacement of the clamped loaded edge. Computational simulation 
indicated a damage initiation load of 4555. 06N (1024 lbs) for the lap joint with epoxy resin. The 
damage initiation mode was by ply longitudinal tensile failures in the 0° ply and the microstress 
damage indexes are activated for the SM1A(+), SM2A(+), SM2B(-), SM3A(+), SM3B(-), SM12A, 
SM12B, SMI 3 A, and SM23A microfailure modes. The locations of damage initiation were at the lap 
joint area. After damage initiation, longitudinal tensile failures spread to the interface between lap joint 
and composite plate as the applied load reached 7584.35N (1705 lbs). Then, longitudinal tensile 
failures continuously spread at the lap joint with the load increasing to 7691 N (1729 lbs). The 
microstress damage indexes were also activated for the SM1A(+), SM2A(+), SM2B(-), SM3A(+), 
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SM3B(-), SMI 2 A, SM12B, SM13A, and SM23A microfailure modes. Thereafter, longitudinal tensile 
failures developed at the lap joint and the interface between lap joint and the composite plate as the 
load increased. The lap joint continued to fracture rapidly and broke under the 44.096kN (9913 lbs) 
loading. Computational simulation indicated that adhesive yielding resulted in failure of the lap joint. 

For the lap joint with graphite/epoxy prepreg tape as its adhesive, the damage initiation load 
was 36.45kN (8193 lbs) by ply longitudinal tensile failures in the 0° ply. The damage initiation took 
place at the edges of the lap joint adjacent to the composite plate. The microstress damage indices were 
activated for SM1A(+), SM2A(+), SM2B(-), SM3A(+), SM3B(-), SM12A, SM12B, SMI 3 A, and 
SM23A microfailure modes. After damage initiation, longitudinal tensile failures spread to plies of the 
specimen as the applied load reached 42.52kN (9558 lbs). Subsequently, damage growth was limited 
mainly to the first two plies of interface between the edge of lap joint and the composite plate until the 
load was increased to 44.53kN (10010 lbs). At 44.75kN (10060 lbs), longitudinal tensile failures 
spread to plies around the edge of lap joint. The microstress damage indices were activated for 
SMI A(+), SM2A(+), SM2B(-), SM3A(+), SM3B(-), SM12A, SM12B, SM13A, and SM23A 
microfailure modes. With increasing load, fracture continued to develop at the edges of the lap joint 
and the specimen broke under the 52.49kN (11800 lbs) loading. Similar to the simulation for the lap 
joint bonded with epoxy resin, adhesive yielding caused failure of bond joint. 

Figure 1.5 shows the relation between displacement and loading for lap joint with two different 
adhesives. It indicates the bonded strength of lap joint with graphite/epoxy prepreg tape is 
approximately 1.19 times higher than that of the epoxy bonded joint, a little lower than experimental 
results by Quattlebaum (1997). This is due to some difference in the material properties in simulation 
and experiment. On the other hand, the apparent linear relationship between the load and displacement 
depicted in Figure 1.5 is not able to reflect the presence of internal damage initiation and growth 
processes. 

For epoxy and graphite/epoxy prepreg tape bonded specimens. Figure 1 .6a and Figure 1 .6b 
compare the change of microstress energy component SM1A(+) with increasing displacement. It 
shows that the peak of SMI A(+) for graphite/epoxy prepreg tape is quite later and higher than that of 
epoxy. This means cracking in region A for epoxy bonded joint is earlier than prepreg bonded joint. 
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The difference in the damage energy amplitudes indicate that the prepreg bonded joint has a much 
higher energy of SMI A(+) damage. 

Figure 1.6c and Figure 1.6d show the microstress energy component SM2A(+) as a function of 
displacement. From the plots, it is obvious that prepreg specimen can withstand higher transverse 
tensile stress in large displacement than epoxy specimen. It reflects that there is greater resistance 
against crack initiation in matrix for prepreg specimen. 

Figure 1.6e and Figure 1.6f show the relation of microstress energy component SM2B(-) versus 
displacement. It is observed that the stress required to cause transverse debonding at the fiber-matrix 
interface for prepreg specimen is larger than that for epoxy specimen. The bond strength for epoxy 
specimen is not as strong as that for prepreg specimen. This is mainly due to the fiber bridging effects 
that impede crack propagation in the prepreg-bonded specimen. 


Figures 1.6g and Figure 1.6h plot microstress energy component SM13A versus displacement. 
In the graph, the shear stress (out of plane) in region A for prepreg bonded joint is higher than that for 
epoxy bonded joint. It shows that there exist bridging fibers in prepreg specimen. 

Recently Nondestructive Evaluation (NDE) methods have been received more attention. Raju 
(1996) studied the in-situ process of C-C composites with Acoustic-Ultrasonic method. Quattlebaum 
(1997) presented the acoustic activity data obtained from epoxy resin and prepreg specimens with the 
same configuration as the computational simulation. Direct comparison of the acoustic emission 
signals with computed microstress level damage energies was not possible because the computed 
damage energy values had much fewer points than the experimentally measured acoustic emission 
data. The apparent reason for the mismatch of refinement between experimental and computational 
values is that microstress level damage tracking is not sufficient to capture the acoustic emission data 
from NDE testing. During the proposed continuation of this project a nanoscale level damage energy 
tracking will be implemented using progressive refinement of the composite subregions via new 
software that have been implemented in the ICAN-Java code. ICAN-Java is a new generation of 
composite mechanics code that enables the telescopic multi-scale subdivision of composites in both the 
thickness and the 3-D spatial sense. Integration of ICAN-Java refinement capability in the 
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CODSTRAN progressive damage and fracture code will be necessary to render progressive 
decomposition of the computationally simulated damage energies to NDE signal precision. For 
correlation of current microstress level damage energy tracking with Quattlebaum’s test data, we used 
a curve fitting software (TableCurve2D v5 2000) to find the fitted curve that envelops acoustic signals. 
Figure 1.7a shows the envelope curve for epoxy resin bonded specimen and Figure 1.7b plots the 
envelope curve for prepreg bonded specimen. Experimental acoustic emission levels are marked with x 
as points through which curves are drawn in Figures 1.7a and 1.7b. The area under the envelope 
represents the total damage energy detected by the acoustic emission during the period monitored. 
Thus, we can correlate the microstress damage energy through computational simulation with the 
relative total damage energies represented by the envelope areas. The envelope area in Figure 1 .7a is 
1.22986e-3 and that in Figure 1.7b is 1.30145e-3. The ratio of the energies from Figure 1.7b to that of 
Figure 1.7a is 1.058. Comparatively, the ratio of simulation results for damage energy shown in Figure 
1.8 is 1.168, which is 10 percent higher than the ratio computed from the NDE test results. 


1.5 CONCLUSIONS 

The failure patterns of the investigated fiber composite specimens and the available computational 
simulation results are found that: 

1. Microstress level damage tracking is able to evaluate the bond strength and monitor adhesive 
yielding. 

2. Microstress level computational simulation provides the details of damage initiation, growth, 
and subsequent fracture in composites. It represents a new approach for investigating damage 
mechanisms of composites. 

3. Computational simulation shows a good correlation with AU signals in the cumulative sense. 

4. More refined nanoscale damage energy tracking via progressive scale decomposition 
implemented in the ICAN-Java code will be necessary to enable direct comparison with NDE 
AU signals. 

5. The demonstrated procedure is flexible and applicable to all types of constituent materials, 
structural geometry, and loading. Hybrid composites, as well as laminated, stitched, woven, and 
braided composites can be simulated. 
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6. Computational simulation, with the use of established composite mechanics and finite element 
models, can be used to predict the influence of microstresses, as well as loading and material 
properties on the durability of composite structures. 
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TABLE 1.1: AS-4 Fiber Properties 


Number of fibers per end = 10000 
Fiber diameter = 0.00508 mm (0.200E 3 in) 

Fiber Density = 4.04E 7 Kg/m 3 (0.063 lb/in 3 ) 

Longitudinal normal modulus = 226.84 GPa (3.29E +7 psi) 

Transverse normal modulus = 13.72 GPa (0.1 99E 7 psi) 

Poisson's ratio ( v n ) = 0.200 
Poisson's ratio ( v 23 ) = 0.250 
Shear modulus (G 12 ) = 13.79 GPa (0.20E +7 psi) 

Shear modulus (G 23 ) = 6.89 GPa (0.1 0E +7 psi) ^ q 

Longitudinal thermal expansion coefficient = 1.0E' 6 /°C (-0.55E 6 /°F) 
Transverse thermal expansion coefficient = 1.0E' 6 /°C (0.56E" 6 /°F) 
Longitudinal heat conductivity = 301 kJ-m/hr/m 2 /°C (4.03 BTU-in/hr/in 2 /°F) 
Transverse heat conductivity = 30.1 kJ-m/hr/nT/°C (0.403 BTU-in/hr/in /°F) 
Heat capacity = 0.712 kJ/kg/°C (0.17 BTU/lb/°F) 

Tensile strength = 3.72 GPa (540 ksi) 

Compressive strength = 3.35GPa (486 ksi) 


TABLE 1.2: Epoxy Matrix Properties 

Matrix density = 3.30E" 7 Kg/m 3 (0.0443 lb/in 3 ) 

Normal modulus = 3.45 GPa (500 ksi) 

Poisson's ratio = 0.35 

Coefficient of thermal expansion = 0.77E' 4 / °C (0.428E' 4 / °F) 

Heat conductivity = 0.648 kJ-m/hr/m 2 / °C (0.868E 2 BTU-in/hr/in 2 / °F) 
Heat capacity = 1 .047 KJ/Kg/ °C (0.25 BTU/lb/ °F) 

Tensile strength = 68.99 MPa (10.0 ksi) 

Compressive strength = 241.59 MPa (35.0 ksi) 

Shear strength = 89.7 MPa (13.0 ksi) 

Allowable tensile strain = 0.02 
Allowable compressive strain = 0.05 
Allowable shear strain = 0.045 
Allowable torsional strain = 0.045 

Void conductivity = 16.8 J-m/hr/m 2 / °C (0.225 BTU-in/hr/in 2 / °F) 
Glass transition temperature = 2 1 6 °C (420°F) 
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FIGURE 1.1: DEFINITIONS OF REGIONS FOR PLY MICROSTRESS CALCULATIONS 



FIGURE 1.2A: UNIT CELL SQUARE ARRAY CONCEPTS OF COMPOSITE 
MICROMECHANICS 
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FIGURE 1 ,2B: HORIZONTAL SLICING FIGURE 1 .2C: VERTICAL SLICING 
FIGURE 1.2: MULTI-SLICE UNIT CELL SUBDIVISIONS 



FIGURE 1.3: TESTING SETUP FOR ADHESIVE BONDED SPECIMEN 



FIGURE 1.4: FINITE ELEMENT MODEL OF COMPOSITE SPECIMEN 
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Loading vs Displacement 



FIGURE 1.5: DISPLACEMENT WITH LOADING FOR COUPON 
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FIGURE 1.6A: DISPLACEMENT WITH SMI A(+) 
FOR EPOXY RESIN BONDED JOINTS 


FIGURE 1.6B: DISPLACEMENT WITH 
SM1A(+) FOR GRAPHITE/EPOXY 
PREPREG BONDED JOINTS 
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FIGURE 1.6C: DISPLACEMENT WITH SM2A(+) 
FOR EPOXY RESIN BONDED JOINTS 
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FIGURE 1.6D: DISPLACEMENT WITH 
SM2A(+) FOR GRAPHITE/EPOXY 
PREPREG BONDED JOINTS 
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FIGURE 1.6E: DISPLACEMENT WITH SM2B(-) 
FOR EPOXY RESIN BONDED JOINTS 


SM2B(-) Fiber-Matrix Interaction Region 



FIGURE 1.6F: DISPLACEMENT WITH 
SM2B(-) FOR GRAPHITE/EPOXY 
PREPREG BONDED JOINTS 




FIGURE 1.6G: DISPLACEMENT WITH SM13A 
FOR EPOXY RESIN BONDED JOINTS 
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FIGURE 1.6H: DISPLACEMENT WITH 
SM13A FOR GRAPHITE/EPOXY 
PREPREG BONDED JOINTS 


FIGURE 1.6: DISPLACEMENT WITH MICROSTRESS LEVEL DAMAGE ENERGY 
COMPONENTS 
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FIGURE 1.7A: THE ENERGY ENVELOPE OF ACOUSTIC EMISSION SIGNALS FOR EPOXY 
RESIN BONDED SPECIMEN 



FIGURE 1.7B: THE ENERGY ENVELOPE OF ACOUSTIC EMISSION SIGNALS FOR PREPREG 
BONDED SPECIMEN 
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FIGURE 1.8: DAMAGE WITH DAMAGE ENERGY FOR COUPON 
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Chapter II. Energy Release Rates and Fracture Planes 


Energy release rate G is usually suggested as a parameter for characterizing material toughness. 
Computational prediction of energy release rate is based on composite mechanics with micro-stress 
level damage assessment, finite element structural analysis and damage progression tracking modules. 
In this chapter, mode I interlaminar and intralaminar energy release rates G, c of composites are 
examined by computational simulation. Results show that computational simulation has a good 
predictive capability for monitoring damage progression in composites. Computational simulation 
enables assessment of the damage initiation and propagation loads. Computational simulation can be 
used prior to testing. Through simulation, sensitive parameters affecting critical values of the energy 
release rates are identified, which significantly enhance the accuracy and productivity of experiments. 
Simulation results are compared with test data. 

Keywords: composites, composite materials, energy release rates, notched beam specimen, 

computational simulation, progressive damage 

2.1 Introduction 

Critical components of a structure are required to remain safe and be able to function under loading 
after experiencing some damage. The cause of damage may be an accident, defect, or unexpected 
overloading. Damage tolerance of a structure is quantified by the residual strength, that is the 
additional load carrying ability after damage. Design considerations with regard to the durability of 
metallic and fiber composite structures require an a priori evaluation of damage initiation and fracture 
propagation mechanisms under expected loading and service environments. Concerns for safety and 
survivability of critical components require a quantification of the structural fracture resistance under 
loading. For a rational design process it is necessary to quantify the structural damage tolerance for a 
candidate design. The ability of designing composites with numerous possible fiber orientation 
patterns, choices of constituent material combinations, ply drops and hybridizations, render a large 
number of possible design parameters that may be varied for an optimal design. The structural fracture 


17 



process of a fiber composite depends on many parameters such as laminate configuration, fiber volume 
ratio, constituent stiffness/strength/hygrothermal parameters, stiffening system, and the fabrication 
process. Recent developments in computational simulation technology implemented in the GRC 
CODSTRAN (Composite Durability STRuctural ANalysis) code have made it possible to evaluate the 
details of progressive damage and fracture in metallic and composite structures. Computational 
simulation enables assessment of the damage initiation and propagation loads. The influence of local 
defects or flaws and effects of the fabrication process in terms of residual stresses are also taken into 
account. 

Computational simulation of 2-D and 3-D fiber composite structures have been successful in 
quantifying the load-displacement relationships, fracture paths, and ultimate strength of composite 
structures. Computational simulation perfonns a complete evaluation of laminated composite fracture 
via assessment of ply and subply level damage/fracture processes. The evaluation of composite test 
response can be made much more productive and informative via computational evaluation of energy 
release rates. Computational simulation can be used prior to testing of a metallic or fiber composite 
structure for evaluation of changes in the local and global critical stress intensity factors and strain 
energy release rates. Progressive damage mechanisms, damage locations/modes, and sensitive 
parameters affecting failure can be identified prior to testing, significantly enhancing the accuracy and 
productivity of an experimental program. 

An important feature of computational simulation is the assessment of damage stability or damage 
tolerance of a structure under loading. At any stage of damage progression, if there is a high level of 
structural resistance to damage progression under the service loading, the structure is stable with 
regard to fracture. The corresponding state of structural damage is referred to as stable damage. On 
the other hand, if damage progression does not encounter significant structural resistance, it 
corresponds to an unstable damage state. Unstable damage progression is characterized by very large 
increases in the amount of damage due to small increases in loading. Whereas during stable damage 
progression the amount of increase in damage is consistent with the increase in loading. 

Internal damage in composites is often initiated as cracking due to normal stresses transverse to 
fiber orientation. At the presence of stress concentrations or defects, initial damage may also include 
fiber fracture. Further degradation is in the form of additional fiber fractures that usually lead to 
structural fracture. Because of the numerous possibilities with material combinations, composite 
geometry, fiber orientations, and loading conditions, it is essential to have an effective computational 
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capability to predict the behavior of composite structures for any loading, geometry, composite 
material combinations, and boundary conditions. The predictions of damage initiation, growth, 
accumulation, and propagation to fracture are important in evaluating the load carrying capacity and 
reliability of composite structures. Quantification of the structural fracture resistance is also required 
to evaluate the durability/life of metallic and composite structures. 

Compared with homogeneous materials, damage initiation and progression characteristics of fiber 
composites are much more complicated. One of the most important problems is interface cracking that 
is known as delamination. Since composite structures often contain some pre-existing or induced flaws 
in matrix and fibers after fabrication of composites, delaminations frequently occur at lower stresses. 
For preventing the degradation in the strength and damage tolerance of a structure due to 
delaminations, it is necessary to quantify the effect of delaminations. 

In general, there are two fracture modes induced by delaminations; interlaminar and intralaminar 
fracture, that are observed in composites. Experimental research [1,2] has been extensively conducted 
on the interlaminar fracture mode. Correspondingly, a few works [3,4] attempted to investigate the 
behavior of intralaminar fracture mode. However, a number of design parameters such as fiber 
orientation patterns, choices of constituent material combinations, ply drops and hybridization, result 
in complex design options for composite structures in which intralaminar transverse tensile fractures 
usually precede more complex fracture modes. Thus, it is difficult for experiments to assess damage 
propagation in different situations. 

Recent developments in computational simulation technology [5,6,7] have made it possible to 
evaluate the details of progressive damage and fracture in composites with the use of established 
material modeling and finite element models, and considering the influence of local defects, through- 
the-thickness cracks and residual stresses. In a computational simulation, damage evolution quantifier 
such as the damage volume, exhausted damage energy, and the damage energy release rate (DERR) are 
used to quantify the structural damage tolerance at different stages of degradation. Low DERR levels 
usually indicate that degradation takes place with minor resistance by the structure. Structural 
resistance to damage propagation is often dependent on structural geometry and boundary conditions 
as well as the applied loading and the state of stress. 

Some simulations [8,9] have been successful in predicting damage tolerance and failure load of 
composite structure by considering ply stresses and the corresponding stress limits for matrix crack 
growth. Thus, sensitive parameters affecting fracture toughness can be identified through simulation, 
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significantly enhancing the accuracy and productivity of experiments. Energy release rate G and 
critical values of the stress intensity factors K are usually suggested as two types of parameters for 
characterizing material toughness. In this chapter, mode I interlaminar and intralaminar energy release 
rates G IC of composites are examined by computational simulation. Results show that computational 
simulation has a good predictive capability for monitoring damage progression in composites. 

Present research models unidirectional Graphite/Epoxy notched beam specimens using the 
CODSTRAN computational code. The effects of the orientation of the fibers with reference to the 
specimen notch direction are investigated with respect to their influences on damage and fracture 
progression characteristics. 

2.2 Methodology 

Computational simulation is implemented via the integration of three modules: (1) composite 
mechanics, (2) finite element analysis, and (3) damage progression tracking. The composite mechanics 
module [5] is designed to analyze fiber composite structures with an updated composite mechanics 
theory. Its main function is to calculate ply and composite properties of laminates from the fiber and 
matrix constituent characteristics and the composite lay-up. Prior to each finite element analysis, the 
ICAN module utilizes a resident data bank that contains the typical fiber and matrix constituent 
properties, computes the composite properties and synthesizes the laminate generalized force- 
displacement relations according to the composite lay-up. Additionally, ICAN determines the 
composite structural response and ply stresses from the FEM analysis results. In simulation, the 
composite mechanics module is called before and after each finite element analysis. The finite element 
analysis module is capable of linear and nonlinear static and dynamic analysis. Four-node anisotropic 
thick shell elements are usually used to model laminated composites [6]. The finite element analysis 
module accepts laminate properties from the composite mechanics module and performs the structural 
analysis at each load increment. After structural analysis, the computed generalized node stress 
resultants and deformations are provided to the composite mechanics module. The composite 
mechanics module computes the developed ply stresses for each ply and checks for ply failure modes 
at each node. Failure criteria applied to detect ply failures are based on the maximum stress criterion 
and modified distortion energy (MDE) for combined stress effects [5]. 
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The overall evaluation of composite structural durability is carried out in the damage progression 
module [9] that keeps track of composite degradation for the entire structure. The damage progression 
module relies on the composite mechanics module for composite micromechanics, macromechanics 
and laminate analysis, and calls the finite element analysis module for global structural analysis. If 
excessive damage is detected, the incremental loads are reduced and the analysis is restarted from the 
previous equilibrium stage. Otherwise, if the increment of loads is acceptable, another finite element 
analysis is performed but the constitutive properties and the finite element mesh are updated to account 
for the damage and deformations from the last simulation. Simulation is stopped when global structural 
fracture is predicted. 


2.3 Microstress Level Damage Tracking 


The matrix in orthotropic composite plies is divided into two parts: regions A and B. Region A 
represents the area in which stress concentrations induced by the interaction of matrix and fiber do not 
create any effect in matrix. Region B represents the interaction zone between fiber and matrix. 
Considering the behavior of longitudinal stress <j iu , transverse stress cr ;22 , in-plane shear stress cr in , 

out-of plane shear stress cr m , temperature gradient AT ; , and moisture M /5 Murthy and Chamis [5] 
present the complete set of equations for evaluating ply microstresses in regions A and B. For example, 
ply microstresses due to cr m are given by: 

^mll = / -^711 )°MI 

CT /1 1 = (^/l 1 / ^71 1 It I 

=(■'.- y ni lEJE m bn, 

(») (b> _„(»> -o-Wl - l_ rr 

~ G mil ~ °/22 “ 33 ; a /ll 


in which cr mll is the matrix longitudinal stress, a fu is the fiber longitudinal stress, <rf, A 22 is the matrix 


transverse stress in region A, and <j^] 2 is the matrix transverse stress in region B. If the ply is subject 

to combined stresses, its microstresses are obtained by simply superimposing results of all 
corresponding stress components. Ply transverse fractures usually begin in region B due to the elevated 
stress levels from stress concentration. Microstress level damage tracking is able to quantify the type of 
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damage in the matrix by comparison of microstresses with constituent stress limits. A microstress 
damage index is defined as a binary number with 14 bits in the damage progression module. The 14-bit 
number corresponds to the following components: (1) SM1A(+) longitudinal stress in region A, 
positive; (2) SMIA(-) longitudinal stress in region A, negative; (3) SM2A(+) transverse stress in region 
A, positive; (4) SM2A(-) transverse stress in region A, negative; (5) SM2B(+) transverse stress in 
region B, positive; (6) SM2B(-) transverse stress in region B, negative; (7) SM3A(+) normal stress in 
region A, positive; (8) SM3A(-) normal stress in region A, negative; (9) SM3B(+) normal stress in 
region B, positive; (10) SM3B(-) normal stress in region B, negative; (1 1) SM12A in-plane shear stress 
in region A; (12) SM12B in-plane shear stress in region B; (13) SM13A out-of-plane shear stress in 
region A; (14) SM23A out-of-plane shear stress in region A. When the binary bit corresponding to 
SM2B(+) is set equal to 1, it indicates that region B is fractured with the transverse failure mode. In a 
subsequent stage transverse cracking will spread to region A. 

2.4 Simulation of Braided Composite Notched Beam Specimens 

Two fiber reinforced bismaleimide prepreg coupons (X5260/G40-800) under three-point bending 
[3] are used in computational simulation. The specimen has a length of L = 27.4mm (1.079 in), a 
width of W = 6.7mm (0.264 in) and a thickness of B = 6.7mm (0.264 in). The fiber volume ratio is 
V { = 0.60 and the void volume ratio is V r = 0.01 . Specimens are labeled as “T-type” and “IF-type” as 

shown in Figures la and lb, respectively. For the IF -type specimen, the crack will grow along the 
fiber direction, whereas the crack will extend orthogonal to the fiber direction in the T -type specimen. 
Specimens are assumed to be dry with zero moisture content. The fiber and matrix properties used for 
computational simulation are given in Table 2.1 and Table 2.2, and the configuration for specimens is 
shown in Figure 2.1a and Figure 2.1b. The initial crack length for T -type specimens is assumed as 
2.96mm (0.117 in) and 3.08mm (0.121 in) for W -type specimen. 

For the IF -type specimen, there are 1275 nodes and 1200 Mindlin type rectangular elements. Its 
finite element model is shown in Figure 2.2b. The notch is modeled to have zero width in the finite 
element model. To represent the test setup, node at one support is restrained against translation to 
model the fixed boundary and node at the other support is constrained only moveable along the length 
of the beam. Load is applied at the top center node. Computational simulation indicated a damage 
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initiation at the crack tip when loading reached 73.40N (16.5 lbs). The damage initiation mode was by 
ply transverse tensile failures in the 0° ply and the microstress damage indices were activated for the 
SMI A(+), SM2A(+), SM2B(+), SM3A(+), SM3A(-), SM3B(+), SMI 2 A, SM12B, SMI 3 A, and 
SM23A microfailure modes. After damage initiation, crack did not extend immediately. When loading 
was applied at 85.27N (19.17 lbs), crack extended again and transverse tensile failures in the 0° ply 
were still the only damage progression mode. The microstress damage indices were activated for 
SM 1 A(+), SM2A(+), SM2B(+), SM3A(+), SM3A(-), SM3B(+), SM12A, SM12B, SMI 3 A, and 
SM23A microfailure modes. With increased loading, transverse tensile failures continuously spread to 
the nodes surrounding the crack tip. The microstress damage indices were also activated for SM1A(+), 
SM2A(+), SM2B(+), SM3A(+), SM3A(-), SM3B(+), SMI 2 A, SM12B, SMI 3 A, and SM23A 
microfailure modes. Then, crack developed quickly with transverse tensile failures occurring in more 
nodes around the crack tip. The specimen continued to fracture rapidly and broke under the 139.54N 
(31.37 lbs) loading. Computational simulation indicated that transverse tensile failures caused fracture 
of the specimen. 

For the T -type specimen, we use 8-node brick elements for constructing the finite element model. 
There are a total of 1977 nodes and 1200 elements. The finite element model is shown in Figure 2.2a. 
Only center node at one support is restrained against translation to model the fixed boundary and the 
other nodes at the same support are constrained only to be moveable in the longitudinal direction of the 
beam. Correspondingly, center node at the other support is constrained only moveable in the 
longitudinal direction and the remaining nodes at the same support are constrained moveable in the 
longitudinal and transverse directions of the beam. Moreover, nodes along the crack edges are tied by 
duplicate node constraints to enforce uniform displacement. The top center nodes in the transverse 
direction are also considered as duplicate nodes to enforce the uniformly displaced loading. The 
damage initiation load was 37.8 IN (8.5 lbs) by ply transverse tensile failures in the 0° ply. The damage 
initiation took place at the crack tip. The microstress damage indices were activated for SM1A(+), 
SM2A(+), SM2B(+), SM3A(+), SM3A(-), SM3B(+), SM12A, SM12B, SM13A, and SM23A 
microfailure modes. After damage initiation, transverse tensile failures spread to nodes around the 
crack tip. As the loading reached 73.04N (16.5 lbs), the crack tip developed additional damage. The 
microstress damage indices were activated for SM1A(+), SM2A(+), SM2B(+), SM3A(+), SM3A(-), 
SM3B(+), SMI 2 A, SM12B, SMI 3 A, and SM23A microfailure modes. Then, crack grew slowly with 
increasing loading. At 101.69N (22.86 lbs), transverse tensile failures spread to more nodes around the 
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crack tip. When the loading reached 106.14N (23.86 lbs), the specimen broke with fracture 
propagation developing very suddenly. The microstress damage indices were also activated for 
SMI A(+), SM2A(+), SM2B(+), SM3A(+), SM3A(-), SM3B(+), SMI 2 A, SM12B, SMI 3 A, and 
SM23A microfailure modes. Similar to the simulation for the W -type specimens, transverse tensile 
fracture mode caused failure of F -type specimen. 

Figure 2.3 shows the relation between displacement and loading for F-type and IF -type 
specimens. It indicates the failure load of F -type is smaller than that of IF -type. From the simulation 
results, the failure load of IF -type is 139.54N whereas it is 106.14N for F -type. Compared with test 
results [3], in which the failure load for IF -type is MON and MON for F-type, it appears simulation 
has a good prediction for final failure load. 

Damage energy and damage volume are another two important metrics for studying damage 
progression in composites. For F -type and IF -type specimens, the relationship between damage 
energy and damage volume is shown in Figure 2.4. It indicates F-type can stand more damage than 
IF -type specimen. Thus, it can be concluded that IF -type is more brittle than F -type specimen. 

Figure 2.5 shows the relationship between load and damage. For IF -type specimen, the damage 
increases slowly until loading reaches 108.98N (24.5 lbs), then damage increases quickly representing 
the spread of transverse tensile failures at more nodes around crack tip. After reaching 132.56N (29.8 
lbs), the damage increase is relatively small until specimen breaks. Whereas, after loading reaches 
11.12N (2.5 lbs), the damage increases uniformly with increased loading for F -type specimen. The 
plot reflects that IF -type specimen is indeed much more brittle than F -type specimen. 

Figure 2.6 shows the relationship between the critical energy release rate Gic and crack area. It is 
observed that Gic for F -type specimen is larger than that for IF -type specimen. This is mainly due to 
the fiber bridging effects that impede crack propagation in the F -type specimen. From the test [3], 
results also show that Gic for F -type specimen is larger than that for IF -type specimen. 

2.5 Conclusions 

Experimental and computational simulation results were compared to investigate the 
progression at crack in composites. The failure patterns of the investigated fiber composite specimens, 
the critical energy release rate G| C and other available computational simulation results indicate that: 
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1. Computational simulation, with the use of established composite mechanics and finite 
element modules, can be used to predict the failure load and fracture toughness of 
composites. 

2. Microstress level computational simulation provides the details of damage initiation, 
growth, and subsequent fracture in composites. It represents a new approach for 
investigating damage mechanisms of composites. 

3. Computational simulation is able to discern damage evolution characteristics with 
regard to fiber and crack orientations. 

4. Fracture toughness characteristics are identified by computational simulation with 
greater perspective compared to only test results. 

5. The demonstrated procedure is flexible and applicable to all types of constituent 
materials, structural geometry, and loading. Hybrid composites, as well as laminated, 
stitched, woven, and braided composites can be simulated. 
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TABLE 2.1: X5260 Fiber Properties 


Number of fibers per end = 1 0000 
Fiber diameter = 0.00762 mm (0.300E’ 3 in) 

Fiber Density = 1743.79 Kg/m 3 (0.063 lb/in 3 ) 

Longitudinal normal modulus = 199.95 GPa (2.9E +7 psi) 

Transverse normal modulus = 19.995 GPa (2.9E +6 psi) 

Poisson's ratio ( v n ) = 0.300 
Poisson's ratio ( v 2i ) = 0.45 

Shear modulus (G 12 ) = 16.548 GPa (2.40E +6 psi) 

Shear modulus (G 23 ) = 8.136 GPa (1.18E +6 psi) 

Longitudinal thermal expansion coefficient = 1.0E' 6 /°C (-0.55E' 6 /°F) 
Transverse thermal expansion coefficient = 1 .0E" 6 /°C (0.56E" 6 /°F) 
Longitudinal heat conductivity - 301 kJ-m/hr/m 2 /°C (4.03 BTU-in/hr/in 2 /°F) 
Transverse heat conductivity = 30.1 kJ-m/hr/m 2 /°C (0.403 BTU-in/hr/in 2 /°F) 
Heat capacity = 0.712 kJ/kg/°C (0.17 BTU/lb/°F) 

Tensile strength = 2.764 GPa (400.9 ksi) 

Compressive strength = 2.41 9GPa (350.9 ksi) 


TABLE 2.2: G40-800 Matrix Properties 

Matrix density = 1264.94 Kg/m 3 (0.0457 lb/in 3 ) 

Normal modulus = 4.669 GPa (720 ksi) 

Poisson's ratio = 0.30 

Coefficient of thermal expansion = 0.77E" 4 / °C (0.428E' 4 / °F) 

Heat conductivity = 0.648 kJ-m/hr/m 2 /°C (0.868E 2 BTU-in/hr/in 2 / °F) 
Heat capacity = 1 .047 KJ/Kg/ °C (0.25 BTU/lb/ °F) 

Tensile strength = 73.706 MPa (10.69 ksi) 

Compressive strength = 55.78 MPa (8.09 ksi) 

Shear strength = 31.16 MPa (4.52 ksi) 

Allowable tensile strain = 0.02 
Allowable compressive strain = 0.05 
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Allowable shear strain = 0.04 
Allowable torsional strain = 0.04 

Void conductivity = 16.8 J-m/hr/m 2 /°C (0.225 BTU-in/hr/in 2 / °F) 
Glass transition temperature = 216 °C (420°F) 
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FIGURE 2.2B: FINITE ELEMENT MODEL (SHELL ELEMENT) FOR W - TYPE SPECIMEN 
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FIGURE 2.3: DISPLACEMENT WITH LOADING FOR IF-TYPE AND T-TYPE COUPONS 
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FIGURE 2.4: DAMAGE ENERGY WITH DAMAGE VOLUME FOR W - TYPE AND 7-TYPE 
COUPONS 
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FIGURE 2.5: LOADING WITH DAMAGE VOLUME FOR IT-TYPE AND T-TYPE COUPONS 



FIGURE 2.6: FRACTURE TOUGHNESS WITH CRACK AREA FOR IT-TYPE AND T-TYPE 
COUPONS 
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Chapter III. Computational Prediction of Effective Elastic Constants 
in a Cross-Ply Laminate under Uniaxial Loading 


The reliability and durability of composite space structures is of critical importance to assure safe 
operation in diverse environments. The main design problem is matrix microcracking at low 
temperatures that adversely affects the durability and increases the permeability of a composite 
structure. Fatigue damage accumulation due to mechanical and thermal cycling causes the growth of 
microcracks into delaminations that cause the failure of interfaces between composite layers. It is 
generally understood that thin plies microcrack less readily than thick plies, and that the mechanical 
strain required to start microcracking depends upon the stiffness of the adjacent "constraint" plies. 
Since microcracking is a function of the difference between adjacent ply angles, a stacking sequence 
chosen for a composite laminate may experience microcracking much more severely than suggested 
from laboratory tests of quasi-isotropic layups. These observations raise concern that the presently 
used methods are not adequate for the design of space structures that will operate reliably at cryogenic 
temperatures. New design methods are needed to take into account the effects of (1) adjacent 
constraint ply properties and fiber angles, (2) ply thickness, and (3) overall ply layup on microcracking 
due to thermomechanical cycling. This report attempts to initiate development of new composite 
mechanics and attendant damage progression/durability codes to take into account the effects of 
adjacent ply stiffness on microcrack initiation and growth. The implementation of this objective 
requires accounting for the stress concentration factors imposed by constraint plies as they affect the 
microcrack growth and evolution. Primary microcracks are typically due to transverse tensile stresses 
and occur as mode 1 dilatation cracks. The new methods account for the effect of ply thickness on 
microcrack formation and durability. The constraints imposed by adjacent plies are shown to be 
significant. The method can be used to model individual fiber-matrix interaction zones and monitor 
changes at sub-constituent levels by subzoning of each sub-slice. The approach can be used to 
evaluate the effects of adjacent ply stress concentrations on microcracking and delamination potentials. 
Thermomechanical cyclic loading efects on microcracking may also be assessed by incorporating the 
constraint ply effects. Initially, the improved method is developed for a cross-ply laminate. It is 
expected that the method can be extended to general angle-plied laminates as well as braided/woven 
composites via the previously demonstrated macromechanics method based on coordinate 
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transformations. It is hoped that the generalization of the method can be carried out in continuation of 
research. 

3.1 Effect of 0° Constraint Plies on Stiffness Degradation in Cross-Ply Laminates 

When cross-ply laminates are subjected to static loading in the 0°ply direction, the first failure 
generally occurs in the transverse plies. Matrix cracking is the main damage mode. Although matrix 
cracking seldom causes the final failure for laminates, it still can greatly impair the performance of 
laminates. Some tests have shown that the stiffness of laminates is reduced because of matrix cracking. 
Moreover, with matrix cracks growing, other damage modes such as delamination and microcracking 
(microcracks extend to plies adjacent to a ply that suffered matrix cracking) will appear. These damage 
mechanisms will lead to laminate failure. Thus, it is necessary to predict the effect of stiffness 
reduction in transverse plies. 

Many analytical models have been developed to evaluate stiffness response of the cracked 
laminate. Berthelot (1997) calculated stress redistribution in cracked cross-ply laminates with assumed 
longitudinal displacement fields. Abdelrahman and Nayfeh (1999) constructed micromechanical 
continuum mixture 2D and 3D models to study the stress redistribution and residual stiffness in 
orthogonally cracked laminates. McCartney (2000) predicted stress distribution in general symmetric 
laminates with uniform cracks by assuming simple displacement patterns that did not satisfy certain 
boundary conditions at the transverse crack-constraint ply interface. 

Pagano and Soni (1983) proposed a global/local model based on plate theory to investigate 
elastic moduli in cracked laminates. Whitney (2000, 2001) determined effective elastic constants of bi- 
directional laminates and angle-ply laminates containing transverse cracks with Pagano’s global/local 
model. Chattopadhyay et al. (1994, 2001) used higher-order plate theory to model transverse matrix 
cracking and local delamination. 

Leblond et al. (1996) developed 2D and 3D numerical models to calculate stiffness reduction in 
cracked cross-ply laminates. Using fracture mechanics and finite element method, Joffe et al (1999, 
2001) analyzed stiffness response in symmetric and balanced laminates due to transverse cracking. 
Whitcomb (2001) employed quasi-3D and 3D finite element methods to evaluate material properties of 
cracked laminates. 
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The variational method based on the principle of minimum complementary energy was first 
applied by Hashin (1985,1987) to study stiffness response of cracked cross-ply laminates under tensile 
loading. In his model, stress components only depend on x-axis direction and are constant across the 
ply thickness. Naim et al. (1989, 1992) used the variational approach to determine residual stress state 
and predict progressive damage in cross-ply laminates. Varna and Berglund (1991) assumed non- 
uniform stress distribution across the ply thickness with variational methods. Praveen and Reddy 
(1998) applied Reddy’s layerwise theory and variational method to analyze stress transfer, stiffness 
reduction and crack opening profiles in cross-ply laminates. Anderssen et al. (1998) also investigated 
stiffness degradation in cross-ply laminates with energy methods, but they assumed displacement 
fields including crack front shape. 

In Anderssen’s model, shear stain was assumed to be zero across 0° plies. Thus, it is obvious 
that the assumption is not realistic and shear stress does not satisfy the continuity condition at the 
interface between 90° and 0° ply. In our research, we improve Anderssen’s model with assuming more 
refined displacement fields. The stress components based on new displacement fields satisfy boundary 
and continuity conditions. The geometry of laminate subjected to tensile loading is illustrated in Figure 
3.1a. The unit cell of the cracked laminate between two adjacent cracks shown in Figure 3.1b is 
adopted for predicting effective elastic constants in the cracked laminate. 
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FIGURE 3.1 A: SCHEMATIC OF A [0°/90°] s LAMINATE WITH TRANSVERSE CRACKS UNDER 
UNIAXIAL LOADING 
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FIGURE 3. IB: ONE-QUARTER OF UNIT DAMAGED CELL WITH TRANSVERSE CRACKING 


3.2 Theoretical Model 

3.2.1 Lamina 1(90°) 

The displacement functions for lamina 1 in Figure 3.1b can be assumed as: 

u ( ?{x,z) = /, {x)+<p{z)f 2 (x) 

U [ ' ] {x,z)= f 2 (x)-y (3-l) 

where f(x),f 2 (x),f} (x) and <fi(z) are unknown functions. 


The boundary conditions for displacement fields are: 


M ( t l) (0,z)= 0 

(3 -2a) 

w i' ) (x.0)= /i(x) 

(3 -2b) 

M^(x,0)= 0 

(3 -2c) 

From (3-2b), we have 


^(o) - o 

(3-3) 

The stress and strain fields are given as: 
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g <o = Hat 


e® = 


dx 

du^ 

dz 


/,) du ^ du ['* 
v' 1 ' = — — + — — 
dz dx 

cr^ = cj'J ■ + c!? • 

r 0) =c 0) . y 0) 

where the elastic stiffness coefficients C{ } are given by 

1 

S 1 1 

E. 


E. 


z 


c .l = 


22 


^ 11*^22 *^12 


5 n .S 2 2 5| 2 


C 66 = 


Substituting (3-1) into (3-4), we find that 


(3-4) 


(3 -4a) 
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4° = // W+^W-Zz W 


4° = 


fA x ) 


=^(z) , /2W + /j'W , f 

*1 

M 

[/, (*) + ^( z ) • fi 4 )] + c u • 

^'(z)-/ 2 (x)+/ 3 (x)4 


a (,) =c (l) 

U z C 12 


r <'W> 

L x z C 66 


(3-5) 


Following the symmetry and traction-free boundary conditions, stress fields are satisfied 


r£ ) (x,0)=0 

r£ ) (0,z)=0 


]<? { ' ] (l,z)dz = 0 

0 

)r^{l, Z yiz = 0 

(3-6) 

Substituting (3-5) into (3-2) and (3-6), we know 

/,(0)=0 

(3 -7a) 

/,(«)=<> 

(3 -7b) 

?>'(«)= 0 

(3 -7c) 

//(o)=o 

(3 - 7d) 

4,)'/ ! W+y/ J '(/) = 0 

On the other hand, it is assumed that 

(3-7e) 

M S 1 ) (°,z)=« S’f 

(3-8) 

Thus, 

A (°) = 

where ' s an unknown constant. 

(3-9) 
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3.2.2 Lamina 2(0°) 

The displacement field is assumed 

u ( p (x,z) = f 3 (*) + m* 2 ' ^ - 

h 

where i/ 0 2 J is an unknown constant. 


(3-10) 


The strain-displacement and stress-strain relations are expressed as: 


tt* 2 ' = 


J2) = 


r? = 


du? 

dx 

dz 

5w* 2 ' 5 m* 2 * 


dz 


dx 


cr* 2 ' = c*?>4 2 ' + c,* 2 'f* 2 ' 
cr* 2 ' = c,* 2 '^ 2 ' + c* 2 '^ 2 ' 


The expression of continuity, symmetry and boundary conditions for lamina 1 and 2 are 
u{ ?( x ’ t ,) =u< ?\ x ' t i) (3 - 1 2a) 

u[^(x,l i )=u^(x,l [ ) (3-1 2b) 

*!?(*.<, )=4, !) (*.'i) (3- 12c) 

<x! l| (*,< 1 )=+ ! +.' 1 ) (3-1 2d) 

tf(0,z)=0 (3 - 1 2e) 

‘fr«(/,z)rfz=0 (3-12f) 


r«»(x,*)-0 (3-12g) 

<+>(*,*) = 0 (3 - 12h) 

Thus, we assume there exists ^(/ ,) and the appropriate expression for r(f(x, z) is 


z™(x,z)=Tf(x,t,)- 
that is 

r* 2) (x,z) = (M, ; (x)- 


h{h + 2/,) 4/7 + 2/, 3 _ 2 

? 2 ^ + Z 

/; /, t; 


(3-13) 


h(h + 2/, ) 4/7 + 2/, 3_^ 2 

- - Z H ~ Z 

/; /, /; 


¥x, W = 4'i • k (', ) ■ fi W + A w] 

Applying (3-10) and (3-1 3a) to (3-1 1) and considering (3-12), we have 


(3 -13a) 
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u?(x, z) = /, (*) + <t>{t, )f 2 (*) -{z-t, )fi (x) 

+ a [^ (/ i )-/ 2 ( x )+ / 3 ' ( x ) 


/;(/; + 2f, ) _ (2/? + 0 2 ^ _1 


2 z ‘+ — z 

t; /, '2 


h 2 t. 


c (l) 
A= 66 


(2) 


(3-14) 


00 

Substituting (3-10) and (3-14) into (3-11) and considering (3-12), we get 

z[ 2) = /,' (*) + <t>{^ )fl (*) - {z - 0 )// (*) 


+ f 2 (x)+ f } " 0 


h{jt + 2t, ) _ (2A + /,) _ 2 | 1 z : 

t 2 t 2 t 2 

*2 l 2 1 2 


// 2 t, 


„( 2 ) 

^.( 2 ) _ M 0 z 

f. 


/*(/? + 2/,) 4A + 2/, ^ 3 

2 2 Z ' 2 ^ 

t; t: t; 


(2) 


r£ ) = ^"( 0 / 2 W+/ 3 / (4 

cr! 2) = cf? [// (x) + <p(t, )/ 2 (x) - (z - 1, )/ 3 (x)] + cf 2 2) 
+ c( 2 UL (ti)-/ 2 (x) + / 3 ' (x)l. 


^ + 2t,) _ (2/» + Q _ 2 | 1 

t; “ t\ t\ 


^ = 4 2) U ( r ) + 4{t\ )fi (*) -(z-t 1 )/," (x)] + C 


,( 2 ) W 0 


(2) 


+ C 




h(h + 2t t ) _ (2A + 0.: . 1 

; i Z +—Z 

t ; /, t; 


h% 

t] 


t] 


P ? = 42 [f (t, )/ 2 (x) + // (x)]- 


From (3-1 2d) and (3-12h), we obtain 

c.?/. (*) + 4 2 V(', )/>' (-0 - 4 2) t 2 /j ' (*) + 4 


/?(/? + 2t,) 4/» + 2t, _ 3 

? ’ 73 z + „2 z 

t; /, t. 


( 2 ) 


( 2 ) K, 


= 0 


t 2 

(2) 


[4> -cg'lr/ M+ [4” -4’R -4 0 

^2 M 

Substituting (3- 16a) into (3-15) and (3-1 6b), we know 


(3-15) 

(3 -16a) 
(3 -16b) 
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of* =f 1 /,'M+‘r,<H',)/(M+kv, -‘tf’jr.'W 

+ c"uU(/ l )-/,(x)+/ J "(.r 


*(* + *,) . (2* + ',) -2 , 1 .» 

/, 2 /j t] 


crj 2) = cgUU' (/,)•/,' (x)+fy (a 


h{h + 2/,)_ (2A + /,) _ 2 | 1 .3 


-z -- 


/; 


t; 


h\ 

t 2 z 

h\ 

tl 


(3 -17a) 


+ c { $(h-z)f?(x) 

<>3fi W+4V(<, )//(*)+«!!’ -^-4 3 U"W= 0 

[>? 

A - J 2 ) l C ' 2 J 

S 1 — C 11 (2) 

c 

^22 

, _M ' 

^ 2 “ „( 2 ) 


(3 -17b) 
(3 -17c) 

(3-1 7d ) 

(3-17e) 


' 22 


3.2.3 Strain Energy 

The strain energy per unit damaged cell is written as: 


U = ^ jj J[cr| l) • f[ 1} + cr[ u ■ s ( ' 1 + • y^Jcfe + j[<x| 2) • f | 2) + a ( : 2) ■ s 


+ c' • k: ] \iz 


(2) ^ r (2)' (2)] 


dx 


■— f u T*(t’ z )d z 

h 


P = hw ■ p {) 
h = /| + f 2 

in which w is the thickness of the laminate. 
Substituting (3-4) and (3-11) into (3-18), we obtain: 


(3-18) 


^ o o 

+^J^ik ) -k! 2) ] 2 


+ 24 ) -4 l) -fl ,) +c 1 (, 1 ) 


+2c 1 (2) -4 ) -4 ) +4 2) 


•k’] 2 + 4 • [ri } ] 2 k 

k 2) ] 2 + c i 2) • \ r »] \ Iz ~ p - «, (2) (0 


0 r, 


(3-19) 


Giving the displacement fields a small variation, then variations of displacement fields and strain 
components are defined as: 
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Su^ (a, z) = Sf (.v) + S<p{z) ■ f 2 (a) + </>(z) ■ Sf 2 (a) 
Su { ; } {x, z) = Sf \ (a) + <p{t, )Sf 2 (a) -(z-t t )5fl (a) 

+ A kio^w+cw]- 

Sui' ) (x,z) = Sf 3 (x)- — 


h{h + 2t { )_ ( 2 /, + 0,2 , 1 ,3 

8 , 2 ~ , 2 
‘2 ‘2 


/ 2 2 /, 


<5w| 2) (a, z) = Sf- j (a) + Su ^ ■ 


izh 


Se ( ' ] = 5f; (a) + 5<f>(z)- f 2 (a) + <f>{z)- Sf 2 (a) 
8e? = Sf 1 (a) + <f>{t , Yfi to W to 


&?) = 


+ 4f(oo7 ! M+<y J 'w]' 

&(*) 


A(A + 2/,)_ (2A + /,)_ 2 , 1 „3 

; * ; z i ' 7 Z 

/; /; t; 


hf , 

? 


&< 2 > = 


( z ) • f 2 (a) + 4 >‘ ’ (2) • Sf 2 (a) + Sf , (a) • 


5 yW = x [( f > (t, (a) + Sf ; (a)] • 


h(h + 2t } ) 4h + 2t A __ 3 


z + — z 
tt t: 


(3-20) 


From (3-3), (3-7a), (3-7b), (3-7e) and (3-9), it is obvious that 

■y,(o)=o 

o /,( 0)-0 

^f,( 0)=Ag 


(3-21) 


<50(O)=O 

The variation of energy £/ is given as: 

su = J<fc}{4’ • 4" • *» + 4 ' • k" • *»' + 4" • &*')+4’ -4' • &!" + C»» ' rS' • SrSh 


0 0 

L h 


r r [cf\ ‘ ^ 2) • &i 2) + c, (2) • [f ( t 2) • <fe! 2) + £'! 2) • <fe (2) ]+ c( 2) • ^ 2) • Se { : } \ 

+ Hl + 4»:4>-^ r 

-P-Su^{l) (3-22) 

Substituting (3-14) into (3-16) and integrating it over z-direction, we have 

SU = SU {]) + SU {2) (3-23) 
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where 


l,U 


*r/0) = f '« '■ '■ 

0 + 1 6 (.v) j"^(z)/z + 1 7 (4 jV(z)^(z)r/z - I g (x) {z)5(t>{z}h 

0 0 o 

I , (x) = 4V, • // (x) + 4> • / 2 (*) • J^(z>/z + 4* • / 3 W 

0 

I 2 M = C 66 • / 2 (*) • ( Z )] dZ + C tb • Au ) ■ // M “ 4 ' 


Z JflZ 


I 3 (x) = 4 •// (*)• j’^(r)fe + 4 ) -/ 2 '(jf)- J^ 2 (z)* + c i ( 2 ) ]</>{z)dz 

\ o 


I 4 W = 4 • /, W + Cl 0 ,’ • + 4 • ^7^ • jV(4 /z 


o 


i s W = 4 • ) • h M - 4 • ^ • j^( z V z + 4', • 


//W 


1 o 








(3 -23a) 


I 7 ( X ) = C 22 ' \f 2 ( X )] 

I 8W= C 66 -|/ 2 W] 2 

and 

"I , {*W; (*)+ I l0 (x)Sf 2 (*)+ 1 , , (xfti (*)+ 1 12 {x)5fi (x)+1 13 (xW (*)' 
_+ 1 ,3 (x)Su^ 


su (2) = J 


dx 

(3 -23b) 
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I 


} W= c f?y 


r,-h ■ /,' W + tfh ■ M>, )• /,'(*)+ 4> %■ *% )• fi W- 4> j ■ f " W 

.2 




0 z 


i ,o M = 4l } • k ('. )] 2 • / 2 M + 4 2 J • <t> (*. ) • fi ( x ) 

I , , (x) = cgVj • ^(/, ) ■ // (x) (/,)•// (x) +c\ 2 ] ) t 2 -f-(t ] )-f; (x) 

3 2 * 3 r , It 2 . . . . , , t 2 


+ c): 


“ “ 12 

1 (J J O ^ 

+ vi;v 3 w- r< 'n 

+ *!?■. 


U 


I W = 4> ^ • <*' (0 • /, M+ 4' ^ • />' W 


hi 


2 

( 2 ) ^ 2 


,W=4 i ^'/'W-4>| 

^- • #' ('| )• /i W + «f? ('.)/! M + 4’ y • // W 

- c <;> 41 . /» ( x >+ c w 41 . // ( x ) _ c«> i- • „«> * <•“> 4 

11 15 105 2 

) = c n ■ fi 


vW'/j'W 


- ..( 2 ) 

15 ' J v ' 11 105 ■" ' ' 2 uz h2 12 0z 

' fi M + c. ( 2 2) • )/a M + 4 2) — 

3/ 1 


L 1Z, 

(2) 

+ c,7 — 


■//W 


Yy-Zj^W+^MY"- “to- 


integrating (3-23a) and (3-23b) again by parts and considering (3-21), it yields 
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- Jfr«2 •/. 


SU = {{[-«, • f"{x)-a 2 •/; (a-)-w, -/j'W-fflj •/ 3 (3) (x)]^r i (4^ 

0 

" (a) + a 3 • f 2 (a) - a 4 • // (.v) + a 5 • // (a) - a 7 • / 3 (3 ) (a )]^/ 2 (a)}</a 

• f '( x )+ G) b -/, (3 ) (a )+« 9 •/ 2 / (A)+<y 10 ■ f } i ) ( x ) + co 0 • /,(a)~L / J 
•/ 3 "(a)+W | 1 -/ 3 (4 ) (a) 

- j{/?, + /?, -<j){z)-P 3 -</>' 

'1 

*i /i (/ ) + <*2/2 W + ©1 /j ( 0‘ ' ' 1 

[(a 2 + a 10 )//(/)+(«,+ a, ,) 


■ 


^,w[ 


+ [«1 fi i l ) + a 2/2 (/)+©, / 3 (/)+ <y 2 // 0 ) + ® 3 «o? - ( z ) 

+ (a, + a 10 )/,' (/ )+ (a 4 + a„ )/: (/) + a 6 / 3 (/ ) + (a 7 +a 12 )//(/) 

+ [- ©./," (/) + ^ 7/2 (0 - ©, 0 / 2 " (/)■ + ©s// (0 - ©,,/j ( 3) (Ok (0 

+ K/i' (0) + ©10/2" (°) + co \ 1 /3 (3, (°)]^S + ["3/1 (0 + ©4/2 w + 
ere 

= c°V +c (2 V 

^22*1 ^ ‘-'ll *'2 

'1 Jit 2 

= 42 jV( Z V Z + c[,V 2 • ^(/, )+ C, (2) (/, ) 

0 

= j[f (z)] 2 + c‘i ) • [f (/, )]" 


*8 W Oz 


7 1V/ 

(«7 +« l 2 )/ 3 , '( / )+« 8 M | (2) 

l +©4/ 2 ( / ) +<y 5 M S ) Ko 2) 


’t-, -<f>{t\)+c^ 

0 0 


+4’<, ^’(0+4’^-k (Of +4’^ 

Ji) '1 -0) '1 932* 

) — }<Kz V z — + (’ ll 1 ~~ • f'('i ) 

*10 ^10 ^ 



o 


,z>/z 


*'(«,)- 4’ 3 I ■ <*'(',)- 4 f ■«>(', ) 
■*'(«,)+ 4 -<K«,)+ 4 - ■<*(», )- 4” £-<*(',) 

r, 6r, 


fe+— JV( Z V 2 -c2J >(?,)- ell’ ■ <t > ' ) 

0 ^ 


'’Jr.M 

(3-24) 
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and 
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(3 - 24c) 


L ~ 

Pi =4 ) JU'M ] dx 

o 


Pi =ci2 JUW ] dx 

0 

According to the theorem of minimum potential energy, there is 
SU = 0 (3-25) 

So, for all admissible displacement, there are 

a, • //' (*) + «2 • // W + °>\ ■ fi ( x ) + 0) 2 • /3 (3) W = 0 ( 3 ■ 26 ) 

«2 ■ /,"(*)- a 3 / 2 W + «4 //W-«5 -//W+«7 ■/ ] <3, ( x ) = ° (3-27) 

co, ■ // (x) + <u 6 ■ /, (3) (x) + a 9 • / 2 ' (x) + co l0 ■ f 2 (3) (x) + co 9 ■ f 3 (x) 

- • fi (*) + 1 ■ /i 4) W = 0 (3-28) 

a+a-*M-a-*"W = 0 ( 3 - 29 ) 

a,// (/) + « 2 fi (0 + ®, / 3 W+ * 2 // (0 + -P = 0 (3-30) 

(a 2 +«lo)/l ( / )+(«4 +a.l)/2 / ( / ) + «6/3( / ) + («7 +a|2)//( / )+« 8 “S? “®12^ = 0 ( 3 ‘ 31 ) 

**/," (/) - *> 7 / 2 (0 - *,// (0+ ®,o// (0 + . /f (0 = 0 (3-32) 

(»)+ ©, 0 // (0) + ^,/, W (0) =0 (3 - 33) 

<*> 3 /1 (0 + «>4 fl (0 + ui 0 z ( 3 - 34 ) 


Solving (3-29), we obtain 

^(z) = ^ sinh(^z) + cosh(^z) - C , 

where 



From (3-3) and (3-7c), there are 

A, = 0 

B , =C, 

Thus, 

<^(z)= CJcosh^zJ-l] 

From (3-1 7c), we know 


(3-35) 

(3 -35a) 
(3 -35b) 

(3 - 36) 
(3-37) 
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/," W = cj : M+ cA<, )■ /,' W + (J> w 




^3 = 


<r 4 = 


c (2 V 

c 12 *2 


,(0 


4 t 

C 12 *1*2 


After substituting (3-38) into (3-26), (3-27) and (3-28), we have 

k J\ (*)+ k 2 fi { x )+ K Jl W = ° 

" {x)-aj 2 {x)+ kJ 2 (x)+Kj; ( a -) = 0 

kJ\ ' (a-) + Kj^ ] (x) + kJ[ (x) + k x J^ ] (a) +K S J 1 (x) = 0 
where 

*•,=«,+ cq 2 £>i 
k 2 =a 2 +co 2 £i -0(0 
k 2 = co x + 6) 2 £ 4 
k 4 =a 2 + a 7 £ 3 
k 5 = a 4 +a 2 C 3 ■ ) 


(3-38) 

(3-39) 

(3-40) 

(3-41) 


K b =a 7 £ 4 -a 5 

k 7 = fi>, +G) U C£ 4 -co s C 3 

*■* =®6 + ^nC 3 


K 9 — ^^11^3^4 ‘^(^l)'®8^3 ’^(0 


*10 = ® 10 +®ll4W(0 


*11 — <i) 9 + ,<1, 4 - 4 

(3 -41a) 

Then, there are two new equations 


ri/.'W + ^/jW+^/a" W = ° 


Ya/\ (a)- r 5 /i (4) (*)+ yJi (a)- y 7 /2 <4) W = 0 

(3-42) 


in which 


V 1 = ^6 

r 2 = ^ 3 «3 

^3 = K 2 K t - * 3*5 

Ya = K i K U ~ K 3 K 1 
Ys = ^ 3 ** 8 

^6 = ^2^11 -^3^9 

r 7 = k 3 k 10 


(3 - 42a) 


From (3-42), we obtain 

Wj} 4) { x )+V / if} 2) ( x ) + V / if2( x ) = 0 ( 3 ‘ 43 ) 

where 
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v'i =r,rs ~rji 

Vi = 7 2 7 5 +7J b -7^4 

¥2 = r 2 r 4 ( 3 ' 43a ) 

Solving (3-43) and considering (3-6) and (3-12) as well as symmetry, we have: 

f 2 {x)- d x sinh(ax)+ J 2 sinh(y 8 r) (3-44) 

for 

^<0 

V'l 

\j/\ - 4^,y / 3 >0 (3 - 44a) 

where d\ and J 2 are constants and 


a = 


P = 


-if/ 2 -‘jyff- 4^, 


^3 


2 ^, 


-V / 2 +W 2 


2y/ x 


>0 


>0 


Substituting (3-44) into (3-42) and considering (3-7), we obtain 

/, (jc) = d 2 x - — + ^ 3 , a t/, -sinh(ca)- ,V: + ^ ■ </, sinh(/&) 

7\P~ 

where ^3 is a constant. 


(3 - 44b) 


(3-45) 


Then, we can know J 3 and /H*) from (3-7), (3-9) and (3-41) 


d 3 = 


+ 




y,a r , 




7\K 7 ^7 

73 


1 


7 y 
3 


X 2 + hP_ + 7dS±lL + 7i K *P K iP K 


7 \P 7y 7^i 


7\ K i 

2 


K 


K 


7 J 




, (i) = +w» +W* -rw -r i*..« [ cosh(ar )_ ,] 

r,K„a 


Mi + r Jl2 £ + J -,* ,/? 1 + r , K,p zL k 1 £_ - r ,*,./? 4 ^ [ cosh ( jSr ) _ ,] + „ o > 


7\ K wP 


(3-46) 


With Newton-Raphson method, we can solve a and (3 from (3-35a), (3-44a) and (3-44b). Then, 
substituting (3-44), (3-45) and (3-46) into (3-30), (3-31), (3-32), (3-33) and (3-34), it is obvious that 
(3-33) is satisfied automatically. Thus, solving equations (3-30), (3-31), (3-32) and (3-34), we obtain 

parameters d\, di, u\'] and u{, 2 J. 
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3.3 Effective Moduli 


All average strain and stress components are calculated with: 
£ (l) = — \dz\g ( ' ) dx 

V 0 0 

g {2) =-^T \dz \g [2) dx 

^ 2 f, 0 

Thus, the effective elastic modulus E x is: 





3.4 Conclusions 


(3-47) 


(3 - 48) 


A new analytical method with spatially consistent continuity conditions has been developed for the 
evaluation of stiffness degradation of cross-ply composites with transverse matrix cracks. The 
developed method is able to quantify the correct residual stiffness of transversely cracked plies, taking 
into account the effects of adjacent constraint ply properties. The development of this method is an 
essential step to improve accuracy of progressive fracture simulations of composite space structures. 
The method can be extended to general angle plied composites via macromechanics and coordinate 
transformations. 
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Chapter IV. Calibration and Validation of Software program ICAN/JAVA 


4.1 Objective 

The objective of this chapter is to present the comparison of results using the old ICAN (4) version and 
the new ICAN/JAVA (8) version. The research also presents the comparison of results using the Multi 
Factor Interaction model not available in the old ICAN. 

4.2 Background 

The most cost effective way to analyze/design composite structures is through the use of computer 
codes. Over the last two decades the research in composite micromechanics and macromechanics, 
which includes the effects of temperature and moisture, has resulted in the development of several 
computer codes for composite mechanics and structural analysis. 

The need for a multilevel analysis to design structural components made of fiber composites is 
necessary for better design. The Multilayered analysis consists of (1) Micromechanical theories for the 
thermoelastic properties and the stress level limit of the single ply as function of constituent material 
properties and the particular fabrication process, (2) the combined stress-strength criteria for the single 
ply, and (3) Multilayered composite structural response and analysis where the interply layered effects 
are taken into account. This code is identified as MFCA (1). 

Intraply hybrid composites are a logical sequel to conventional and interply hybrid composites. 
Another code INHYD (2) integrated the mechanical behavior of Intraply hybrids with those for 
hygrothermal effects for predicting hygral, thermal, and mechanical properties of Intraply hybrid 
composites for their design. 

ICAN is the synergistic combination of the micro mechanical design of INHYD and the laminate 
analysis of MFCA with several significant enhancements. It was primarily designed to analyze the 
hygro-thermo-electro-mechanical properties and response of fiber- or particulate-reinforced, resin- 
matrix or metal-matrix layered composites, given the local membrane loads and bending moments. 
Three types of layers are recognized by the program: (1) the standard composite system that consists 
entirely of a primary composite made of one type of fiber and matrix; (2) the intraply hybrid composite 
system that consists of a primary composite and a secondary composite arranged in a prescribed 
manner within a layer (For purposes of identification, the primary composite in the hybrid is the one 
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that constitutes the largest volume ratio.); and (3) the interply layer that consists of the matrix. In 
addition, ICAN recognizes moisture, temperature, and electrical gradients through the thickness. 
However, within each ply (or slice) the temperature, moisture, or voltage is assumed to be constant. 
During the last six years, since its intial release ICAN has undergone many changes to fix the bugs that 
have surfaced and also improvement to its micromechanical equations. ICAN/JAVA is the result of 
these enhancements which had made it more manageable and user-friendly. 

4.3 Features and Enhancements of ICAN/JAVA 

(1) The output can now be tailored to specific needs by choosing the appropriate options. 

(2) Several modules have been added to perform durability/fatigue type analyses for thermal as well as 
mechanical cyclic loads. The code can currently assess degradation due to mechanical and thermal 
cyclic loads with or without a defect. 

(3) The laminate configuration is not restricted to only plies but can be sliced and subsliced for a closer 
look at what goes on in the ply layer. 

(4) The laminate configuration can also be used to incorporate any defects in the layup. 

(5) The thermal loads, hygral loads and electrical loads can now be input as constant, linear, parabolic, 
hyperbolic or user defined across the ply-layup. 

(6) The damping analysis has been incorporated in the new version. 

(7) Details regarding the reaction of the metal and matrix composite can also be input in the new 
version. 

(8) Details regarding the impact of a hard spherical projectile crashing into the composite can also be 
input in this version. 

(9) The ability to change empirical constants like the experimental correlation factor for combined 
stress value criteria, and for onset of delamination criteria value is also provided. 


4.4 Theory 

In ICAN/JAVA the initiation of damage is detected using modified distortion energy (MDE) 
failure criteria (33) The MDE failure criteria is a variation of the Tsai-Hill theory (34,35) that states 
that failure is initiated when the following inequality is violated. 


F = 1- 
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V ^Li2s J 


<0 


(4.1) 


Where K[ [2 is the longitudinal-transverse directional elastic properties interaction factor. 
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(4.2) 


given in terms of the lamina elastic constants by 

„ _ 0 ~ YLn]_^_L22 0 ~ ^7,23 )^L] I 

A. i f 

JE LU E L22 (2 + V Ln +V'i|3 )(2 + V L 2 ] + ^i23 ) 
oui and, ol 22 are the lamina stresses along and perpendicular to the fibers and 01 . 12 s is the shear stress 
and F is the failure. 

Slii and Sl 22 , are the lamina strengths along and perpendicular to the fibers and Sli2s is the shear 
strength. If the prevailing stress in the fiber direction is tensile (oui a 0) then the strength S L n is the 
tensile longitudanal strength Slut and if the prevailing stress in the fiber direction is compressive 
(auia 0) then the strength Slii is the compressive longitudanal strength Sluc- The same is applied for 
the transverse stresses. 

Eli i and El 22 are the lamina moduli along and perpendicular to the fibres 
Vli 2 , v L i 3 and Vl 23 are the lamina poisson’s ratios 

Failure is deemed to take place when F=0 and failure modes can be either fiber or matrix failure 
depending upon the dominanat terms in the MDE equation. Fiber failure is assigned when one of the 
following equations is satisfied 
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S Llu/ 

The stress limits Slut, Sluc and Su2s in the MDE failure criteria are calculated by ICAN based on the 
constituent fiber and matrix strength and micromechanics equations. The lamina elastic properties used 
in the directional interaction factor Kli 2 are also computed by ICAN from the fiber and matrix elastic 
properties via micromechanics equations (4). 

The simulation of complex material behavior resulting from the interaction of several factors such as 
temperature, fatigue, time dependence etc has been mainly performed by factor-specific 
representations like creep analysis, fatigue analysis etc. Suppose we assume material behavior to be a 
continum represented by some surface (Referenec 12). Initially a primitive form of MFIM 
representation for describing complex composite behavior in polymer matrix composites (36). This 
was extended to metal matrix composites (37) and is continuing. 
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A multifactor interaction model is briefly described to represent complex point material behavior in a 
single equation. The model is of product form in order to represent coupled interactions and to be 
computationally effective. The model describes a continuum or surface in space that represents the 
complex material behavior in terms of various factors that affect a specified material behavior as 
shown in equation below (38). 
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where M=property, T=temperature, S=strength, o=stress, N=number of cycles, t=time, to=load 
frequency, E e =erosion effect, C e =corrosion effect, C c =chemical/metallurgical effect 
Subscripts: gw=wet glass temperature, o=reference condition-assumed nominal at ambient 
conditions, f=final condition, M=mechanical load, T=thermal cyclic load, M,n,etc are exponents for 
the mterial that property effect which describe respective behavior paths from the reference to the 
final value. 

ICAN/JAVA has integrated the following factors in its code: thermal, stress, stress rate, temperature 
rate, reaction, mechanical cycle fatigue, thermal cycle fatigue, time at stress, moisture, porosity, 
hygrothermal effect and frequency effects. 

Reference value, final value, increasing exponent, decreasing exponent, and other exponent are entered 
in the material bank of ICAN/JAVA. 

The equation of the multiplicative factor is 

- “| exp onent 

final _ value - factor _ value (4 5) 

final _ value - reference value _ 

where the program chooses the exponent in the computation from the three exponents (increasing 
exponent, decreasing exponent and other exponent) as appropriate. For the hygrothermal effect the 
equation of the multiplicative factor is slightly modified as below (Reference 14): 

“I exp onent 

T -T 

(4.6) 

T -T 

1 gd J o 
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where T 0 is the room temperature with a default value of 70 degrees Fahrenheit, T gd is the dry glass 
transition temperature, entered as the final value, T gw is the wet glass transition temperature. T gw is not 
entered but is computed from T gd and the moisture at the current layer as below: 

= T gd (1 - 0. 1 M f + 0.005 Mf ) (4.7) 

where M f is the moisture fraction, T is the computed temperature value of the current ply/layer/slice 


4.5 Model Description 


4.5.1 Durability Analysis 

The laminate consists of 1 6 plies as shown in Figure 4. 1 with a ply thickness of 0.005 inch. 
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-45 DEGREE PLY 

7 

45 DEGREE PLY 

6 

0 DEGREE PLY 

5 

90 DEGREE PLY 
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-45 DEGREE PLY 

3 

45 DEGREE PLY 

2 

0 DEGREE PLY 

l 


Figure 4. 1 : Composite Ply Lay-Up for durability analysis 

The material is made of T300/IMHS composite the properties of which are shown in Chapter 5.10. 
The fiber volume ratio for the primary material is 0.55 while for the secondary material is 0.54. The 
cure temperature (Tcu) is 160F while the use temperature (Tu) is 350F. The moisture content is 0.5%. 
The laminate is subjected to in-plane membrane loads of 100 pounds per inch and 50 pounds per inch 
in the x and y directions. Two different cases are run and the results compared for durability. In 
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addition to the static loads in the first case the tension-tension cyclic load of 200 pounds per inch is 
input. The cyclic degradation coefficient is taken as 0.01 and the number of cycles is 10 million. 

In the second case the tension-tension and the bending cyclic load of 200 pounds inch per inch is input 
with the same cyclic degradation coefficient and number of cycles as shown in Table 4.2. 


Table 4. 1 :Durability Analysis Loads for the test cases 


Mnemonic 

Test Case 1 

Test Case 2 

Nxx 

100 

100 

Nyy 

50 

50 

CNxx 

200 

0 

CMxx 

0 

200 

Cyclic Degradation 

0.01 

0.01 

Coefficient 



Number of Cycles 

l.E+07 

l.E+07 


This problem is compared with ICAN, ICAN/JAVA (subslicing=0) and ICAN/JAVA (subslicing = 9 
). The BOOLEAN MICRO is set true because the ply properties were not found in the online databank 
for ICAN/JAVA. Hence to compare correctly both the MICRO keywords are set to true in the input 
files. 

4.5.2 MFIM Analysis 

The laminate consists of 4 plies as shown in Figure 4.2 with a thickness of 0.01 inch for zero degree 
plies and 0.005 inch for ninety degree plies. 
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Figure 4.2: Composite ply lay-up for MFIM analysis 

It is made of AS— /IMLS and SGLA/HMHS (Chapter 5.10) composite materials with the first applied 
to the zero degree plies and the second to the ninty degree plies. The primary fiber volume ratio for 
AS— /IMLS and SGLA/HMHS is 0.55 and the secondary fiber volume ratio is 0.57. The cure and use 
temperature is 70F. 

The laminate is subjected to in-plane static load of 1000. pounds per inch in the x-direction and a 
combination of cyclic loads shown in Table 4.1 . 
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Table 4.2: Loads api 

plied for MFIM test 

Mnemonic 

Upper limit of 
cyclic load 

Lower limit of 
cyclic load 

Number of cycles 

Cyclic 

degradation 

coefficient 

(Betal) 

CNXX 

200. 

100. 

100. 

0.1 

CNYY 

-50. 

-100. 

10. 

0.1 

CNXY 

20. 

10. 

100. 

0.2 

CMXX 

10. 

5. 

10. 

0.01 

CMYY 

4. 

2. 

1000. 

0.15 

CMXY 

2. 

1 . 

100. 

0.1 

NXX 

1000 

- 

- 

- 


The material properties are further changed with the help of the MFIM using the following format as 
shown in Table 4.3. 


Computed Factor 

Reference 

Value 

Final 

Value 

Increasing 

Exponent 

Decreasing 

Exponent 

Stress (psi) 



0.5 

0.5 

Stress rate (psi/sec) 



0.5 

0.5 

Temperature rate (F/sec) 

73 

73 

0.34 

0.34 

Mechanical Cycle 

Fatigue (hertz) 



0.5 

0.5 

Time at stress (sec) 



0.5 

0.5 


The material properties that were changed were the Poisson’s ratio, the elastic moduli and the 
strengths. The new composite was named AS-1/IML1 and SGL1/HMHS and the same simulation 
carried out and the results compared. 
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4.6 Results and Discussion 

The durability results are compared for the different loads. The graphs are plotted between the safety 
factors and the ply-layup. A brief explanation of these safety factors is as follows. Ml 1 is a safety 
factor defined as below: 


M\\ = \-R\\ 


where R\ 1 = 

S l wr or ^ LUC 

where R1 1 is defined as the ratio of the stress in 1 1 direction to the 
compression. 

M22 is defined below: 


(4.8) 

ultimate stress in tension or 


M22 = \-R22 

whereRll = — ^ 4 ' 9 ^ 

S L22T OI 'S L22C 

where R22 is defined as the ratio of the stress in 22 direction to the ultimate stress in tension or 
compression. 

M12 is a safety factor defined as below: 

Ml 2 = \-R\2 

a n (4.10) 

where R\ 2 = — — 
c 

°L I2S 

where R12 is defined as the ratio of the stress in 12 direction to the ultimate shear stress. 


4.6.1 Cyclic In Plane Membrane Load 

Cyclic in plane membrane loads were input and the results computed from ICAN (Case 1), 
ICAN/JAVA (sub-slicing=0) (Case 2) and ICAN/JAVA (sub-slicing=9) (Case 3) 

The factor of safety for the transverse, longitudinal and shear directions were calculated and 
enumerated in Tables 4. 4-4. 6 and the variations plotted with respect to ply lay-up as shown in Figures 
4. 3-4.5. 


Table 4.4: Factor of safety in the longitudanal direction for cyclic in plane membrane load 




Mil 

PLY 

ANGLE 

ICAN 

ICAN/JAVA 
(Without sub-slicing) 

ICAN/JAVA 
(With sub-slicing) 


57 
















Figure 4.4: graph of m22 for Cyclic In Plane Membrane Load 

Table 4.6: Factor of safety in the shear direction for cyclic in plane membrane load 
PLY I ANGLE f M12 


















From the figures we observe that plies 4 and 13 are the first to fail in the transverse and longitudinal 
directions, which is consistent with both ICAN/JAVA and ICAN. But in case of shear the ±45 plies are 
the first to fail and the 0 and 90 plies do not take the shear loads at all. This is consistent with the loads 
applied. Also we observe that there is a marked difference between the factors of safety for 
ICAN/JAVA and ICAN. This can be due to the sub-slicing used in ICAN/JAVA, which enables the 
code to predict the stresses in more detail and accuracy. Hence the factors of safety are usually higher 
or equal in ICAN/JAVA than in ICAN. 

4.6.2 Cyclic In Plane Bending Load 

Cyclic in plane bending loads were input and the results computed from ICAN (Case 1), ICAN/JAVA 
(sub-slicing=0) (Case 2) and ICAN/JAVA (sub-slicing=9) (Case 3) 

The factor of safety for the transverse, longitudinal and shear directions were calculated and 
enumerated in Tables 4. 7-4.9 and the variation plotted with respect to ply lay-up as shown in Figures 
4. 6-4.8. 


Table 4.7: Factor of safety in the longitudanal direction for In plane bending load 




M22 




ICAN/JAVA 

ICAN/JAVA 

ICAN/JAVA 

PLY 

ANGLE 

ICAN 

(Without sub- 

(With sub-slicing) 

(With sub-slicing) 




slicing) 

Slice 1 

Slice 125 

1 

0 

-4.33 

-1.594 

-1.766 

-1.423 

2 

45 

-0.15 

0.427 

0.384 

0.470 

3 

-45 

-0.32 

0.341 

0.283 

0.399 

4 

90 

0.032 

0.541 

0.488 

0.594 

5 

0 

-1.49 

-0.221 

-0.392 

-0.049 

6 

45 

0.543 

0.768 

0.726 

0.811 

7 

-45 

0.622 

0.807 

0.749 

0.865 

8 

90 

0.927 

0.965 

0.912 

0.982 

9 

90 

0.811 

0.911 

0.981 

0.840 

10 

-45 

0.760 

0.888 

0.931 

0.844 

11 

45 

0.701 

0.859 

0.891 

0.827 

12 

0 

-0.84 

0.117 

0.246 

-0.011 

13 

90 

-0.39 

0.344 

0.415 

0.273 

14 

-45 

0.055 

0.538 

0.582 

0.495 

15 

45 

0.185 

0.603 

0.635 

0.571 
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Figure 4.7: graph of m22 for Cyclic In Plane Bending Load 
Table 4.9: Factor of safety in the shear direction for In plane bending load 
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10 

-45 

0.675 

0.717 

0.811 

0.623 

11 

45 

0.463 

0.529 

0.623 

0.434 

12 

0 

0.953 

0.958 

0.964 

0.952 

13 

90 

0.940 

0.946 

0.952 

0.939 

14 

-45 

- 0.17 

- 0.037 

0.057 

- 0.131 

15 

45 

- 0.39 

- 0.226 

- 0.132 

- 0.320 

16 

0 

0.900 

0.909 

0.915 

0.903 



Figure 4.8: graph ofm12 for Cyclic In Plane Bending Load 

From the figures we observe that the plies (0/45/-45/90/0) from the bottom and (0/90/-45/45/0) top fail 
in case of application of bending loads in the longitudinal direction, which is expected while the 
transverse direction has the plies (90/-45/45/0) from the top failing. In case of shear the plies (-45/45) 
fail from the ends followed by the (-45/45) plies in the center. Also we observe that in case of all the 
predictions for slice 1 and slice 125 the ICAN/JAVA with no sub-slicing predicts the average of the 
two values, which is consistent throughout. 

4.6.3 MFIM Results 

The MFIM is a new addition to ICAN/JAVA and a durability analysis was carried out to observe the 
effects of MFIM on the results. Figures 4.9-4. 1 1 show the plots for the various factors of safety. 
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—♦—WITHOUT MFIM 
-•-WITH MFIM 


1 (0 DEG) 


2 (90 DEG) 3 (90 DEG) 4 (0 DEG) 

PLY-LAYUP 


Figure 4.9: graph of Ml 1 for MFIM 



1 (0 DEG) 2 (90 DEG) 3 (90 DEG) 4 (0 DEG) 

PLY-LAYUP 


Figure 4.10: graph of m22 for MFIM 





Figure 4.11: graph of m 1 2 for MFIM 

In Figure 5.16, the stresses in the 1 1 directions vary in their values by only 5%. The difference is 
reduced to about 0.5% in Figure 5.17 for the case of stresses in 22 directions. While the shear stresses 
show a difference of 16% in the center of the ply-layup. This shows that MFIM indeed has some effect 
on the results for the same models and loads. An attempt was made to try to bring the two variations as 
close by using different coefficients for MFIM modeling. We can see that the 0 degree plies fail first in 
shear and transverse directions while in case of longitudinal direction the 0 degree ply in the bottom 
fails first. 

4.7 Conclusion 

From the above comparisons we can observe that there are some differences in the results of the 
software ICAN and ICAN/JAVA. This may be due to the addition of new modules like sub-slicing, 
variable loading etc in ICAN/JAVA. The other cause can be some bugs that need to be fixed. We also 
observe that the MFIM has an effect on the results as the material properties are changed to a 
considerable extent. This is important because then we can take into account the change in material 
properties with application of loads over time to get more accurate results. 

4.8 Future Work 

Hence, a further detailed look into the new ICAN/JAVA for its corroboration with the ICAN results is 
required to fix any bugs that might be present in the new software. Also a need to further study in 
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detail the new enhancements of ICAN/JAVA for better results is necessary to realize the full potential 
of the software. 
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Chapter V 

Prediction of First Ply Failure and Fracture in Composite Materials of Different 
Size and Geometry 

5.1 Objective 

The objective of this chapter is to simulate the biaxial failure of composites of different size and 
geometry and compare them with experimental results. The research consists of generation of first-ply 
failure envelopes for combined loading of these laminates on the basis of first-ply failure and laminate 
fracture. 

5.2 Background 

For the last three decades or so, efforts have been made to predict the laminate fracture under uniaxial 
and/or combined loading. Due to the existence of a large number of failure criteria and analysis 
methods no exact solution has been found which can correctly predict the composite laminate fracture. 
But a comparison of the various failure theories with a benchmarked experimental data can be useful 
and instructive for structural designers. 

To that end, Hinton and Soden (27) organized an exercise to compare the current theories of composite 
failure with their experimental results. Among the many theories Gotsis et al. (28) used the 
micromechanics-based theory and progressive fracture. The results were comparable with the 
experimental data but not very accurate for the biaxial failure envelopes. Others like Sun et al. (29) 
used the linear laminate analysis, Rotem used the Rotem failure criteria (30) and Worlfe et al. (31) 
used the strain-energy based criteria to compare with the experimental data. Some of the theories could 
accurately predict the stress/strain curves while others could predict the biaxial failure but no one 
theory could predict the laminate fracture for all the benchmarks and hence only a comparison of many 
theories occurred. An attempt has been made to predict the biaxial failure envelope for composites 
using the CODSTRAN and ICAN in this research and to compare with the data of Hinton et al. (27) 

5.3 Theory 

CODSTRAN is used to simulate multi-layered fiber composite structures under any loading and 
hygrothermal conditions. It permits a simulation of local behavior as well as global structural behavior 
as shown in Figure 5.1 
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Figure 5.1 : Simulation cycle of progressive fracture in composite laminates and 
STRUCTURES VIA CODSTRAN 

The constituent properties are described as functions of environmental and mechanical loading 
conditions as shown in the lower part of the figure. By knowing the developed ply stresses, and ply 
strength we can predict the damage initiation, growth, accumulation, and propagation, which results in 
the constituent material properties being updated at every load increment. For example, if a ply’s 
transverse strength exceeds its allowable strength then the ply is assumed to fail and the matrix 
modulus is updated as negligible. But if a ply longitudinal stress exceeds the allowable strength, then 
the fiber/matrix longitudinal moduli are both updated as negligible. This means that the ply is assumed 
to carry no load and the stresses are then redistributed to the surrounding plies. Once the current 
constituent properties are determined the repeated application of micromechanics, macromechanics, 
and laminate theory are used to assemble the global structural stiffness matrix, which is fed to the finite 
element analysis. 

The nodal stress resultants are from the finite element analysis and are used to decompose the changes 
in the global response (laminate stress and strains) on the local material stress/resistance. The load is 
increased only if there is no further damage due to changes in ply level stresses. Otherwise, only the 
material properties are updated at every increment until equilibrium is reached between the applied 
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load and the local response. Until a specified convergence is reached the structural equilibrium is 
maintained as shown in Figure 5.2. 



Figure 5.2: CODSTRAN load increment 

The final result in terms of load as a function of global displacement is shown in Figure 5.3. Agarwal 
et al. (32) have described the damage progression analysis procedure in detail. 



Figure 5.3:Overall CODSTRAN simulations 
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5.4 Model Description 

Two different specimens were tested under biaxial stress state to generate the failure envelopes. The 
material used was E-glass/MY750 [±55] s (Appendix B) laminate. A square one-inch by one-inch 
coupon was tested under transverse and shear loading as shown in Figure 5.4. The coupon model was 
similar to the one used by Gotsis et al. (28) to represent the stress state prior to the disclosure of the 
cylindrical specimen and it’s boundary conditions, by Hinton and Soden (27). 



Figure 5.4: E-glass/MY750 [± 55] s Coupon 

The results were further compared with a cylindrical laminate subjected to a uniform internal pressure. 
The length of the specimen is 8.4 inches and the tube is clamped at both ends. The inner diameter of 
the tube is 39.4 inches with a thickness of 0.4 inches as shown in Figure 5.5. 
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Figure 5.5: E-glass/MY750 [± 55] s Cylindrical Tube 

End reinforcement with variable thickness made of E-glass/MY750 [90] is applied to the tube. The 
reinforcement is provided to simulate real testing conditions where the specimen is held between the 
ends. The main focus is at the center of the tube for a length of 23.6 inches where the E-glass/MY750 
[ ±55] s laminate is present with no reinforcement. There laminate lay-ups used the base of [±55] s as 
the basic core and the reinforcements were added as 90° lamina each of thickness 0.01578 inch. To 
simulate the variable reinforcement twelve ply-layup were defined and were applied to different 
sections of the cylinder as shown in the Figure 5.6. Hence, the thickness of the cylinder varied from 
0.039 inches at the center to 0.25 at the ends as shown in Figure 5.7. 
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Figure 5.7: Variable Thickness of the cyclinder due to different ply-schedules 
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The cylindrical tube was further tested with the calibrated properties SGE1/MHD1 and the results 
compared with the non-calibrated material properties. These material properties are shown in 
Appendix B. 

5.5 Model Assumptions 

The coupon was fixed in the z direction along the boundary. The first node at the origin was fixed in 
all directions while the x-axis and y-axis were fixed in y and x directions respectively. The thickness of 
the coupon was 0.039 inches, which is the same as the thickness of the center part of the cylinder. 
Axial compression and tension were applied to the coupon along the x and y directions to simulate the 
stress state. A force of 1 kip was used as the initial load increment. 



Figure 5.8: Boundary Conditions for the coupon 

To simulate the boundary conditions of clamped specimen for the cylinder one end of it was fixed in 
all directions while the other end was fixed in y and z directions only as shown in Figure 5.8. This 
allowed for deformation of the specimen in the axial direction. 
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z 


X 



Figure 5.9: E-Glass/MY750 Tube Geometry 

The axial force and internal pressure were applied as loads to simulate the experimental conditions as 
shown in Figure 5.9 and the results plotted. 



Figure 5.10: Direction of force and pressure applied to the cylinder specimen 
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The stress ratio was varied from —1 to 1 and the following procedure used to calculate the required 
internal pressure. 


(7 „ = 


a v = 


EL 

t 

YA 

2Urt 


StressRatio = — - 
cr 


(5.1) 


where, 

p - int ernal _ pressure 
r = radius _of _ cyclinder 
t - thickness of _ cyclinder 

The force applied was 1 kip and the total force on the top of the cylinder was 32 kips. We calculated 
the internal pressure for SR=2/1 as follows. 


2Urtr 


YZi _ 32 

nr 2 n(1.97) 2 


2.62 psi I inch 2 


(5.2) 


where, 

p = pressure 

7, F t = summation _of _ forces 
r = radius _of _ cyclinder = 1 .91 inch 
t = thickness _of _cyclinder = 0.039 inch 

5.6 Material Calibration 

AS CODSTRAN is micromechanics based code, a calibration of the material properties so as to match 
the predicted composite unidirectional properties were undertaken. The following graph shown in 
Figure 5.10 was used to calibrate the properties. 
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Figure 5.11: Biaxial Envelope of 0° GRP Lamina under Combined a x and a Y Stresses for E- 
glass/MY750 Epoxy. 

ICAN was used to carry out the material calibration as it has the ability to calculate the laminate 
properties from the fiber and matrix properties. The model used for the calibration was a four ply 0° 
unidirectional laminate. The existing material properties of SGES/MHDY (Appendix B) were used as 
the basis of our calibration. In-plane loads were applied to check the response to the variation of the 
material properties. 

Table 5.1 shows the comparison of the stresses obtained from simulated and corrected material 
properties to that of experimental data. 


Table 5.1 : Experimental and simulated data for material calibration 


PMEMB 


STRESS 

(Mpa) 

STRESS 

(Mpa) 



EXPERIMENTAL 


SIMULATED 

IN X- 

IN Y- 

SIGMA X 

SIGMA Y 

SIGMA X 

SIGMA Y 

DIRECTION 

DIRECTION 





+200 

0 

+ 1271.86 

0.000 

+ 1290.545 

0.000 

-200 

0 

-830.508 



0.000 

0 

+ 10 

0.000 

+42.373 

0.000 

+66.647 

0 

-10 

0.000 

-140.678 

0.000 

-99.975 
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From Figure 5.1 1 we observe the difference in the biaxial failure envelope for the lamina to be less than 
7% in the sigma x direction and less than 40% in the sigma y direction with respect to the experimental 
data. 


Experimental Simulated 



Figure 5. 1 2:Comparison of Biaxial Envelope of 0° GRP Lamina under Combined ct x and o y 
Stresses for experimental and simulated material properties. 

The error in the y direction was higher because a slight change in properties could increase the error of 
the x direction and this was the closest we could get to match the envelope overall. Tables 5.2 and 5.3 
show the new material properties and Chapter 5.10 has the original properties. Comparisons of the 
fabric and matrix properties are shown in Table 5.4 and 5.5. 


Table 5.2: SGE1 -Fiber Properties 

SGE1 Silenka E- GLASS 1 200 tex FIBER (Soden). 


$ 

$ 

$ 

Number of fibers per end 


Nf 3000 number 
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Filament equivalent diameter df 

0.300E-03 

inches 

Weight density 

Rhof 

0.640E-01 lb/in* *3 

Normal moduli (11) 

Efl 1 

1.073E+07 

psi 

Normal moduli (22) 

Ef22 

1.073E+07 

psi 

Poisson’s ratio (12) 

Nufl2 

0.200E+00 

non-dim 

Poisson’s ratio (23) 

Nuf23 

0.200E+00 

non-dim 

Shear moduli (12) 

Gfl2 

4.467E+06 

psi 

Shear moduli (23) 

Gf23 

4.467E+06 

psi 


Thermal expansion coef. ( 1 1 ) Alfafl 1 0.272E-05 in/in/F 
Thermal expansion coef. (22) Alfaf22 0.272E-05 in/in/F 
Heat conductivity (1 1) Kfl 1 5.208E-02 BTU-in/hr/in 2 /F 

Heat conductivity (22) Kf22 5.208E-02 BTU-in/hr/in 2 /F 

Heat capacity Cf 0.170E+00 BTU/lb/F 

Fiber tensile strength SfT 3.073E+05 psi 
Fiber compressive strength SfC 2.103E+05 psi 


Table 5.3: MHD1 -Matrix Properties 
MHD1 My7 5 0/H Y 9 1 7/D Y063 matrix (Soden) 

$ 

$ 

$ 

Weight density Rhom 0.470E-01 lb/in**3 

Normal modulus Em 0.486E+06 psi 

Poisson"s ratio Num 0.350E+00 non-dim 

Thermal expansion coef. Alfa m 0.322E-04 in/in/F 
Matrix heat conductivity Km 8.681E-03 BTU-in/hr/inVF 

Heat capacity Cm 0.250E+00 BTU/lb/F 

Matrix tensile strength SmT 0.737E+04 psi 

Matrix compressive strength SmC 0.245E+05 psi 
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Matrix shear strength SmS 0.122E+05 psi 

Allowable tensile strain eps mT 0.500E-01 in/in 

Allowable compr. strain eps mC 0.500E-01 in/in 
Allowable shear strain eps mS 0.500E-01 in/in 

Allowable torsional strain eps mTOR 0.500E-01 in/in 
Void heat conductivity kv 0.225E+00 BTU-in/hr/in 2 /F 

Glass transition temperature Tgdr 0.420E+03 F 



5.7 Results And Discussion 

The simulations were done in three parts. Fist the various uniaxial and biaxial loads were applied to the 
coupon and the results for the first-ply failure were converted to the stresses in x and y directions and 
tabulated in Table 5.6. 
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Table 


APPLIED INITIAL FORCE 

RESULTANT 

FAILURE 

FAILURE STRESS (Mpa) 

INCREMENT (kip) 

FORCE (kip) 




IN X- 

IN Y- 

IN X- 

IN Y- 

SIGMA X 

SIGMA Y 

DIRECTION 

DIRECTION 

DIRECTION 

DIRECTION 



+i 

0 

+0.4429 

0.0000 

+078.30 

0.00 

-i 

0 

-0.3907 

0.0000 

-69.07 

0.00 

0 

+1 

0.0000 

+1.4770 

0.00 

+261.10 

0 

-1 

0.0000 

-3.1630 

0.00 

-559.18 

-i 

-1 

-1.9560 

-1.9560 

-345.80 

-345.80 

+i 

+1 

+ 1.6540 

+ 1.6540 

+292.40 

+292.40 

+i 

-1 

+2.2310 

-2.2310 

+394.40 

-394.40 

-i 

+1 

-0.4387 

+0.4387 

-077.55 

+077.55 


The following formula was used to convert the force into stress for the coupon. 


a = 


Lt 


* 6.89415Mpa 


(5.3) 


where, 

L = length of _ coupon = 1 .0 inch 
t = thickness of _ coupon = 0.039 inch 


Secondly the non-calibrated properties of SGES/MHDY laminate were used to simulate the 


biaxial/uniaxial loading and tabulated in Table 5.7. 

Table 5.7 Cylinder Simulation Data: SGES/MHDY Laminate 


INITIAL 

FORCE 

INCREMENT 

IN X- 

DIRECTION 

(kip) 

APPLIED 

INTERNAL 

PRESSURE 

(psi) 

RESULTANT 
FAILURE 
FORCE IN X- 
DIRECTION 
(kip) 

RESULTANT 

FAILURE 

PRESSURE 

(psi) 

SIGMA X 
(Mpa) 

SIGMA Y 
(Mpa) 

+i 

0.0000 

+5.842 

0.0000 

+83.44 

0.00 

-i 

0.0000 

-6.615 

0.0000 

-94.98 

0.00 

+i 

-1.3123 

+3.655 

-0.1499 

+52.20 

-52.20 

-i 

+1.3123 

-4.513 

+0.1851 

-64.46 

+64.46 
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+1 

+ 1.3123 

+7.179 

+0.2944 

+102.54 

+ 102.53 

-1 

-1.3123 

-9.725 

-0.3988 

-138.91 

-138.89 

0 

+ 1.3123 

0.000 

+0.5024 

0.00 

+ 174.97 

0 

-1.3123 

0.000 

-0.4061 

0.00 

-141.43 

+1 

+2.6246 

+ 10.930 

+0.8961 

+156.12 

+312.08 

+1 

+2.2730 

+ 11.830 

+0.8400 

+168.97 

+292.55 

+1 

+4.5931 

+7.280 

+1.0450 

+103.98 

+363.94 


Thirdly the calibrated properties were subjected to the same loads and the results tabulated in Table 5.8 


Table 5.8: Cylinder Simulation Data: SGE1/MHD1 Laminate 


INITIAL 

FORCE 

INCREMENT 

IN X- 

DIRECTION 

(kip) 

APPLIED 

INTERNAL 

PRESSURE 

(psi) 

RESULTANT 
FAILURE 
FORCE IN X- 
DIRECTION 
(kip) 

RESULTANT 

FAILURE 

PRESSURE 

(psi) 

SIGMA X 
(Mpa) 

SIGMA Y 
(Mpa) 

+i 

0.0000 

+4.565 

0.0000 

+65.21 

0.00 

-i 

0.0000 

-8.379 

0.0000 

-119.68 

0.00 

0 

+1.3123 

0.000 

+0.5450 

0.00 

+ 189.80 

0 

-1.3123 

0.000 

-0.3024 

0.00 

-105.32 

+1 

+1.3123 

+5.190 

+0.2129 

+74.13 

+74.14 

+1 

-1.3123 

+2.513 

-0.1030 

+35.89 

-35.87 

-1 

+ 1.3123 

-5.469 

+0.2243 

-78.11 

+78.11 

-1 

-1.3123 

-13.720 

-0.5626 

-195.97 

-195.94 

-1 

-3.2808 

-7.744 

-0.7940 

-110.61 

-276.53 

+1 

+3.2808 

+7.929 

+0.8129 

+113.25 

+283.11 

+1 

+4.2808 

+6.833 

+0.9141 

+97.60 

+318.35 

+1 

+8.2808 

+2.288 

+0.5920 

+32.68 

+206.17 

+1 

+0.2808 

+4.727 

+0.0415 

+67.52 

+ 14.440 
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+1 

+4.2809 

+6.844 

+0.9155 

+97.75 

+318.84 

+1 

+4.2812 

+6.833 

+0.9142 

+97.60 

+318.39 

+1 

+4.3000 

+5.783 

+0.7771 

+82.60 

+270.64 

+1 

+4.2830 

+6.703 

+0.8971 

+95.74 

+312.43 

+1 

+4.2815 

+6.740 

+0.9018 

+96.27 

+314.07 

+1 

+4.2811 

+6.182 

+0.8270 

+88.30 

+288.02 

+1 

+4.2812 

+6.833 

+0.9142 

+97.60 

+318.39 

+1 

+4.2800 

+6.029 

+0.8063 

+86.11 

+280.81 

+1 

+4.2810 

+6.844 

+0.9156 

+97.75 

+318.88 


The following formula was used to convert the force into stress for the cylinder. 

V f 

cr Y =±±^-*6.89475Mpa (5.4) 

v 2U rt 

where, 

r = radius of _ cyclinder = 1 .91 inch 
t = thickness of _cyclinder = 0.039 inch 

The following formula was used to convert the pressure into Mega Pascal. 

cr K =^*6.89475 Mpa (5 5) 

Figure 5.12 and 5.13 show the comparison of the experimental and theoretical results for the biaxial 
failure envelopes of E-glass/MY750 [±55] s laminate, coupon and the SGE1/MHD1 [±55] s laminate 
under combined o x - a y biaxial loading. 
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Sigma X (MPa) 


—♦—Coupon (SGES/MHDY) — Cylinder (SGES/MHDY) Cylinder (SGE1/MHD1) 



Figure 5.13: Biaxial Failure Stresses for [ ±55] s E-Glass/MY750 Laminate using 
CODSTRAN 
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Figure 5.14: Biaxial Failure Stresses for [±55] s E-Glass/MY750 Laminate Obtained 
Experimentally 

The theoretical final strength computed by CODSTRAN under predicted the experimental leakage 
strength between SR=1:1 and SR=2:1. In the compression-compression quadrant CODSTRAN 
predicted that the composite had no damage tolerance or the ultimate fracture to occur immediately 
after first ply failure due to longitudinal compressive stresses. These predictions were approximately 
three times lower than the test results at a stress ratio of approximately —2:1. This may be attributed to 
(a) the high fiber volume fraction of the test specimens used in this quadrant [16], (b) the use of thick 
tubes in the experiments while the theory uses thin flat plates, (c) simulations neglected the residual 
matrix stiffness of the failed plies that may have under predicted the ultimate failure load under 
compressive stresses. The final strength for the coupon as predicted by CODSTRAN is higher than the 
cylindrical specimen because the coupon is not subjected to internal pressure but a combination of o x - 
a y stresses. The additional boundary conditions used alongside the edges of the coupon, were 
responsible for the stiffening of the coupon resulting in higher values for ultimate failure load. There 
were no experiments for the coupon and only the cylinder results were compared with the given 
experimental data. 
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5.8 Conclusion 

A comparison was made between the theoretical predictions and progressive fracture simulations using 
CODSTRAN and it was found that only a limited number of segments of the biaxial failure envelope 
for [±55] s glass/epoxy laminate correlated to the test results. The final strengths were under-predicted 
in a vast majority of segments. This can be attributed to the high fiber volume fraction of the test 
specimens, the use of thick tubes in the experiments while the simulations use a thin plate and the 
simulations neglected the residual matrix stiffness of the failed plies that may have predicted the 
ultimate failure loads under compressive stresses. 

The ability of CODSTRAN to predict first ply failure envelopes is fairly accurate provided the 
structural condition of the specimen is known. 

5.9 Future Work 

A need to simulate the exact structural conditions at the time of experimentation to compare with the 
experimental results is required. The different models like thin and thick tubes need to be taken into 
account in the simulation to get more accurate comparison. Residual matrix stiffness after transverse 
tensile failure needs to be quantified and accounted for in the damage simulations. 
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5.10 Material Properties 

Fiber Properties: 


AS-4 


Number of fibers per end = 10000 
Fiber diameter = 0.00762 mm (0.300E-3 in) 

Fiber Density = 4.04E-7 Kg/m 3 (0.063 lb/in 3 ) 

Longitudinal normal modulus = 227 GPa (32.90E+6 psi) 

Transverse normal modulus = 13.7 GPa (1.99E+6 psi) 

Poisson's ratio (v I2 ) = 0.20 
Poisson's ratio (V 23 ) = 0.25 
Shear modulus (G !2 ) = 13.8 GPa (2.00E+6 psi) 

Shear modulus (G 23 ) = 6.90 GPa (1.00E+6 psi) 

Longitudinal thermal expansion coefficient = -1.0E-6/°C (-0.55E-6 /°F) 
Transverse thermal expansion coefficient = 1.0E-5/°C (0.56E-5 /°F) 
Longitudinal heat conductivity = 43.4 J-m/hr/m 2 /°C (580 BTU-in/hr/uT/°F) 
Transverse heat conductivity = 4.34 J-m/hr/nT/°C (58 BTU-in/hr/inV°F) 
Heat capacity - 0.712 KJ/Kg/°C (0.17 BTU/lb/°F) 

Tensile strength = 3.723 GPa (540 ksi) 

Compressive strength = 3.351 GPa (486 ksi) 

SGES SilenkaE- GLASS 1200 tex FIBER (Soden). 


Number of fibers per end 

Nf 

3000 number 

Filament equivalent diameter df 

0.300E-03 

inches 

Weight density 

Rhof 

0.640E-01 lb/in* *3 

Normal moduli (11) 

Efl 1 

1.073E+07 

psi 

Normal moduli (22) 

Ef22 

1.073E+07 

psi 

Poisson”s ratio (12) 

Nufl2 

0.200E+00 

non-dim 

Poisson"s ratio (23) 

Nuf23 

0.200E+00 

non-dim 

Shear moduli (12) 

Gfl2 

4.467E+06 

psi 

Shear moduli (23) 

Gf23 

4.467E+06 

psi 


Thermal expansion coef. (11) Alfafl 1 0.272E-05 in/in/F 
Thermal expansion coef. (22) Alfaf22 0.272E-05 in/in/F 
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Heat conductivity (11) Kfl 1 

5.208E-02 

BTU-in/hr/in**2/F 

Heat conductivity (22) Kf22 

5.208E-02 

BTU-in/hr/in**2/F 

Heat capacity Cf 0. 

170E+00 BTU/lb/F 

Fiber tensile strength SfT 

3.1 18E+05 

psi 

Fiber compressive strength SfC 

2.103E+05 psi 


T300 GRAPHITE FIBER. 



Number of fibers per end 

Nf 

3000 number 

Filament equivalent diameter df 

0.300E-03 

inches 

Weight density 

Rhof 

0.640E-01 lb/in**3 

Normal moduli (11) 

Efl 1 

0.320E+08 

psi 

Normal moduli (22) 

Ef22 

0.200E+07 

psi 

Poisson"s ratio (12) 

Nufl2 

0.200E+00 

non-dim 

Poisson"s ratio (23) 

Nuf23 

0.250E+00 

non-dim 

Shear moduli (12) 

Gfl2 

0.130E+07 

psi 

Shear moduli (23) 

Gf23 

0.700E+06 

psi 


Thermal expansion coef. ( 1 1) Alfafl 1 -.550E-06 in/in/F 
Thermal expansion coef. (22) Alfaf22 0.560E-05 in/in/F 
Heat conductivity (11) Kfl 1 0.403E+01 BTU-in/hr/in**2/F 

Heat conductivity (22) Kf22 0.403E+00 BTU-in/hr/in**2/F 

Heat capacity Cf 0.170E+00 BTU/lb/F 

Fiber tensile strength SfT 0.350E+06 psi 
Fiber compressive strength SfC 0.300E+06 psi 


AS- GRAPHITE FIBER. 


Number of fibers per end 

Nf 

10000 

number 

Filament equivalent diameter df 

0.300E-03 

inches 

Weight density 

Rhof 

0.630E-01 lb/in* *3 

Normal moduli (11) 

Efl 1 

0.310E+08 

psi 

Normal moduli (22) 

Ef22 

0.200E+07 

psi 

Poisson M s ratio (12) 

Nufl2 

0.200E+00 

non-dim 

Poisson"s ratio (23) 

Nuf23 

0.250E+00 

non-dim 

Shear moduli (12) 

Gfl2 

0.200E+07 

psi 

Shear moduli (23) 

Gf23 

0.100E+07 

psi 
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Thermal expansion coef. (11) Alfafl 1 -.550E-06 in/in/F 
Thermal expansion coef. (22) Alfaf22 0.560E-05 in/in/F 
Heat conductivity ( 1 1 ) Kfl 1 0.403E+01 BTU-in/hr/in**2/F 

Heat conductivity (22) Kf22 0.403E+00 BTU-in/hr/in**2/F 

Heat capacity Cf 0.170E+00 BTU/lb/F 

Fiber tensile strength SfT 0.400E+06 psi 
Fiber compressive strength SfC 0.400E+06 psi 

SGLA S- GLASS FIBER. 

Number of fibers per end Nf 204 number 


Filament equivalent diameter df 

0.360E-03 

inches 

Weight density 

Rhof 

0.900E-01 lb/in**3 

Normal moduli (11) 

Efl 1 

0.124E+08 

psi 

Normal moduli (22) 

Ef22 

0.124E+08 

psi 

Poisson"s ratio (12) 

Nufl2 

0.200E+00 

non-dim 

Poisson"s ratio (23) 

Nuf23 

0.200E+00 

non-dim 

Shear moduli (12) 

Gfl2 

0.517E+07 

psi 

Shear moduli (23) 

Gf23 

0.517E+07 

psi 


Thermal expansion coef (1 1) Alfafl 1 0.280E-05 in/in/F 
Thermal expansion coef (22) Alfaf22 0.280E-05 in/in/F 
Heat conductivity ( 1 1 ) Kfl 1 5.208E-02 BTU-in/hr/in**2/F 

Heat conductivity (22) Kf22 5.208E-02 BTU-in/hr/in**2/F 

Heat capacity Cf 0.170E+00 BTU/lb/F 

Fiber tensile strength SfT 0.360E+06 psi 
Fiber compressive strength SfC 0.300E+06 psi 
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Matrix Properties: 


5250 HMHS 


Matrix density = 3.40E-7 Kg/m 3 (0.0457 lb/in 3 ) 

Normal modulus = 4.62 GPa (671 ksi) 

Poisson’s ratio = 0.705 

Coefficient of thermal expansion = 0.518E-4/°C (0.288E-4 /°F) 

Heat conductivity = 0.649E-3 J-m/hr/rrT/°C (0.862E-8 BTU-in/hr/inV°F) 
Heat capacity = 1 .047 KJ/Kg/°C (0.25 BTU/lb/°F) 

Tensile strength = 75.1 MPa (10.9 ksi) 

Compressive strength = 283 MPa (41.0 ksi) 

Shear strength = 138 MPa (20.0 ksi) 

Allowable tensile strain = 0.02 
Allowable compressive strain = 0.05 
Allowable shear strain = 0.04 
Allowable torsional strain = 0.04 

Void conductivity = 16.8 J-m/hr/irf/°C (0.225 BTU-in/hr/in"/°F) 

Glass transition temperature = 216°C (420°F) 

MHDY My750/HY917/DY063 matrix (Soden) 


Weight density Rhom 0.470E-01 lb/in**3 

Normal modulus Em 0.486E+06 psi 

Poisson”s ratio Num 0.350E+00 non-dim 

Thermal expansion coef. Alfa m 0.322E-04 in/in/F 

Matrix heat conductivity Km 8.681E-03 BTU-in/hr/in**2/F 

Heat capacity Cm 0.250E+00 BTU/lb/F 

Matrix tensile strength SmT 0.149E+05 psi 

Matrix compressive strength SmC 0.174E+05 psi 

Matrix shear strength SmS 0.100E+05 psi 

Allowable tensile strain eps mT 0.500E-01 in/in 

Allowable compr. strain eps mC 0.500E-01 in/in 

Allowable shear strain eps mS 0.500E-01 in/in 

Allowable torsional strain eps mTOR 0.500E-01 in/in 

Void heat conductivity kv 0.225E+00 BTU-in/hr/in**2/F 

Glass transition temperature Tgdr 0.420E+03 F 
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IMLS INTERMEDIATE MODULUS LOW STRENGTH MATRIX. 


Weight density Rhom 0.460E-01 lb/in**3 

Normal modulus Em 0.500E+06 psi 

Poisson"s ratio Num 0.410E+00 non-dim 

Thermal expansion coef. Alfa m 0.570E-04 in/in/F 

Matrix heat conductivity Km 8.681E-03 BTU-in/hr/in**2/F 

Heat capacity Cm 0.250E+00 BTU/lb/F 

Matrix tensile strength SmT 0.700E+04 psi 

Matrix compressive strength SmC 0.210E+05 psi 

Matrix shear strength SmS 0.700E+04 psi 

Allowable tensile strain eps mT 0.140E-01 in/in 

Allowable compr. strain eps mC 0.420E-01 in/in 

Allowable shear strain eps mS 0.320E-01 in/in 

Allowable torsional strain eps mTOR 0.320E-01 in/in 

Void heat conductivity kv 0.225E+00 BTU-in/hr/in**2/F 

Glass transition temperature Tgdr 0.420E+03 F 

IMHS INTERMEDIATE MODULUS HIGH STRENGTH MATRIX. 

Weight density Rhom 0.440E-01 lb/in* *3 

Normal modulus Em 0.500E+06 psi 

Poisson"s ratio Num 0.350E+00 non-dim 

Thermal expansion coef. Alfa m 0.360E-04 in/in/F 

Matrix heat conductivity Km 8.681E-03 BTU-in/hr/in**2/F 

Heat capacity Cm 0.250E+00 BTU/lb/F 

Matrix tensile strength SmT 0.150E+05 psi 

Matrix compressive strength SmC 0.350E+05 psi 

Matrix shear strength SmS 0.130E+05 psi 

Allowable tensile strain eps mT 0.200E-01 in/in 

Allowable compr. strain eps mC 0.500E-01 in/in 

Allowable shear strain eps mS 0.350E-01 in/in 

Allowable torsional strain eps mTOR 0.350E-01 in/in 
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Void heat conductivity kv 0.225E+00 BTU-in/hr/in**2/F 

Glass transition temperature Tgdr 0.420E+03 F 

HMHS HIGH MODULUS HIGH STRENGTH MATRIX. 

$ 

$ 

$ 

Weight density Rhom 0.450E-01 lb/in* *3 

Normal modulus Em 0.750E+06 psi 

Poisson"s ratio Num 0.350E+00 non-dim 

Thermal expansion coef. Alfa m 0.400E-04 in/in/F 

Matrix heat conductivity Km 8.681E-03 BTU-in/hr/in**2/F 

Heat capacity Cm 0.250E+00 BTU/lb/F 

Matrix tensile strength SmT 0.200E+05 psi 

Matrix compressive strength SmC 0.500E+05 psi 

Matrix shear strength SmS 0.150E+05 psi 

Allowable tensile strain eps mT 0.200E-01 in/in 

Allowable compr. strain eps mC 0.500E-01 in/in 

Allowable shear strain eps mS 0.400E-01 in/in 

Allowable torsional strain eps mTOR 0.400E-01 in/in 

Void heat conductivity kv 0.225E+00 BTU-in/hr/in**2/F 

Glass transition temperature Tgdr 0.420E+03 F 
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